None Ch2_exe
In [138]:
print(R.version)
file.path(R.home("bin"), "R")
               _                           
platform       x86_64-apple-darwin13.4.0   
arch           x86_64                      
os             darwin13.4.0                
system         x86_64, darwin13.4.0        
status                                     
major          4                           
minor          3.1                         
year           2023                        
month          06                          
day            16                          
svn rev        84548                       
language       R                           
version.string R version 4.3.1 (2023-06-16)
nickname       Beagle Scouts               
'/Users/karlzhang/miniforge3/envs/JTR/lib/R/bin/R'
In [139]:
library(IRdisplay) # display
library(fUnitRoots) # adfTest
library(TSA) # eacf
In [109]:
# detach("package:AFTSCode", unload = TRUE) # First detach the package
# unloadNamespace("AFTSCode")
In [140]:
library(AFTSCode)
In [ ]:
# plot_auto_arima_forecast_fig <- function(
#     da_ts, eotr, h, npts, frequency,
#     xreg=NULL, main=NULL, xlab=NULL, ylab=NULL, ylim = NULL, ...
# ) {
#     par(bg = "white")
#     # arima model
#     if (is.null(seasonal)) {seasonal = list(order = c(0L, 0L, 0L), period = NA)} # default value
#     tr_da_ts <- ts(da_ts[1:eotr], frequency = frequency, start = start(da_ts))
#     if (is.null(xreg)) {
#         tr_xreg <- NULL
#         fc_xreg <- NULL
#     } else {
#         stopifnot("xreg should be of type matrix"=(is.matrix(xreg)))
#         stopifnot("length(da_ts)!=dim(xreg)[1]"=(length(da_ts)==dim(xreg)[1]))
#         tr_xreg <- xreg[1:eotr]
#         fc_xreg <- xreg[(eotr+1):dim(xreg)[1]]
#     }
#     ts_fm <- auto.arima(tr_da_ts, xreg = tr_xreg, ...)
#     print(ts_fm$nobs)
#     # Forecast
#     ts_fm$x <- tr_da_ts # https://stackoverflow.com/a/42464130/4307919
#     ts_fc_res <- forecast(ts_fm, h = h, xreg = fc_xreg)
#     # Plot forecast
#     if (is.null(npts)) {
#         npts <- eotr
#     }
#     xmin <- time(da_ts)[eotr] - npts / frequency
#     xmax <- time(da_ts)[eotr] + (max(h, length(da_ts) - eotr) + 1) / frequency
#     cat(xmin, ";", xmax)
#     # # Label 1: Actual Observation Line
#     plot(ts_fc_res, xlim = c(xmin, xmax), ylim = ylim, main = main, xlab = xlab, ylab = ylab)
#     # Plot forecast mean (prepend the last observed data in the training dataset)
#     dummy_1st_fmean_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$mean)), frequency = frequency, start = end(tr_da_ts))
#     # # Label -: NULL
#     lines(dummy_1st_fmean_ts)
#     # # Label 2: Forecast Mean
#     points(dummy_1st_fmean_ts, pch = 1)
#     # Plot confidence interval (95%)
#     dummy_1st_flower_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$lower[, 2])), frequency = frequency, start = end(tr_da_ts))
#     dummy_1st_fupper_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$upper[, 2])), frequency = frequency, start = end(tr_da_ts))
#     # # Label 3: Forecast 95% Lower Bound
#     lines(dummy_1st_flower_ts, lty = 2)
#     # # Label 4: Forecast 95% Upper Bound
#     lines(dummy_1st_fupper_ts, lty = 2)
#     # Plot original data
#     orig_plot_ts <- ts(da_ts[(eotr - npts + 1):length(da_ts)], frequency = frequency, start = time(da_ts)[eotr] - (npts - 1) / frequency)
#     # # Label -: NULL
#     lines(orig_plot_ts)
#     # # Label 5: Actual Observation Points
#     points(orig_plot_ts, pch = 19)
#     legend(
#         "topleft", 
#         legend = c(
#             "Actual Obs Line", NULL, "Forecast Mean", "Forecast 95% Lower Bound",
#             "Forecast 95% Upper Bound", NULL, "Actual Obs"
#         ), 
#         # col = c("black", "red", "blue"), 
#         lty = c(1, NA, 2, 2, NA),
#         pch = c(NA, 1, NA, NA, 19)
#     )
#     ts_fc_res
# }

2-1

  • From Sec. 2.5.4
    • $E[R_{101}|F_{100}] = c_0 - \theta_1 a_{100} = 0.2*0.01 = 0.002$
    • $Var[R_{101}|F_{100}] = \sigma_{a}^2 = 0.025^2 = 0.000625$
    • $E[R_{102}|F_{100}] = c_0 = 0$
    • $Var[R_{102}|F_{100}] = (1+\theta_1^2)\sigma_{a}^2 = (1+0.2^2)*0.025^2 = 0.00065$
  • From Sec. 2.5.1
    • $\rho_1 = -(-0.2)/(1+0.2^2) = 0.192307692307692$
    • $\rho_2 = 0$

2-2

  • From Sec. 2.4.1
    • $E[r_t]=\frac{\phi_0}{1-\phi_1}=0.01/(1-0.2)=0.0125$
    • $Var[r_t]=\frac{\sigma_a^2}{1-\phi_1^2}=0.02/(1-0.2^2)=0.0208333333333333$
  • From Sec. 2.4.2
    • $\rho_0=1$
    • $\rho_1=\phi_1\rho_0=0.2$
    • $\rho_2=\phi_1\rho_1=0.04$
  • From Sec. 2.4.4
    • $E[r_{101}|F_{100}]=\phi_0+\phi_1*r_{100}=0.01+0.2*(-0.01)=0.008$
    • $Var[r_{101}|F_{100}]=\sigma_a^2=0.02$
    • $E[r_{102}|F_{100}]=\phi_0+\phi_1E[r_{101}|F_{100}]=0.01+0.2*0.008=0.0116$
    • $Var[r_{102}|F_{100}]=Var[\phi_1 r_{101}+a_{102}|F_{100}]=Var[\phi_1 (\phi_0+\phi_1 r_{100}+a_{101})+a_{102}|F_{100}]=(1+\phi_1^2)\sigma_a^2=(1+0.2^2)*0.02=0.0208$

2-3

In [10]:
da = read.table("../AFTS_sol/data/m-unrate.txt", header = T)
da[1:5,]
A data.frame: 5 × 4
YearMonDayRate
<int><int><int><dbl>
11948113.4
21948213.8
31948314.0
41948413.9
51948513.5
In [15]:
unem_rate = da$Rate
lg_unem_rate = log(unem_rate)
unem_rate_ts = ts(unem_rate, frequency = 12, start = c(1948,1))
lg_unem_rate_ts = ts(lg_unem_rate, frequency = 12, start = c(1948,1))
In [143]:
plot_time_fig(unem_rate_ts, main = "Unemployment Rate", xlab = "Year")
pacf(unem_rate)
No description has been provided for this image
No description has been provided for this image
In [22]:
plot_time_fig(lg_unem_rate_ts, main = "log(Unemployment Rate)", xlab = "Year")
pacf(lg_unem_rate)
No description has been provided for this image
No description has been provided for this image

Unit-root test

  • We can reject the null-hypothesis from ADF test. In other words, the process doesn't have unit-root.
  • Seems like we cannot reject he null hypothesis when we use large enough lags in the ADF testing.
  • Also from the ACF plots of the unemployment rate before and after taking difference, we can see that we need to take difference to cancel the strong serial correlation.
  • From the ACF plot of the unemployment rate after taking difference, we see some damping sine and cosine wave, indicating that there could be some business cycle (pp42).
In [34]:
# "lags=12" is chosen from PACF plot
adf_test_res = adfTest(unem_rate_ts, lags = 12, type = c("ct"))
adf_test_res
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -2.9715
  P VALUE:
    0.1671 

Description:
 Tue Feb 27 23:20:00 2024 by user: 
In [31]:
lag.max = 60
plot_acf(da = unem_rate, lag.max = lag.max, main = "Unemployment Rate")
pacf(unem_rate, lag.max = lag.max)
No description has been provided for this image
No description has been provided for this image
In [35]:
diff_unem_rate = diff(unem_rate)
diff_unem_rate_ts = ts(diff_unem_rate, frequency = 12, start = c(1948,2))
In [36]:
plot_acf(da = diff_unem_rate, lag.max = lag.max, main = "diff(Unemployment Rate)")
pacf(diff_unem_rate, lag.max = lag.max)
No description has been provided for this image
No description has been provided for this image
In [144]:
plot_acf(da = diff(unem_rate, 12), lag.max = lag.max, main = "diff(Unemployment Rate, 12)")
pacf(diff(unem_rate, 12), lag.max = lag.max)
No description has been provided for this image
No description has been provided for this image

Order determination

  • Try to use ARIMA(4, 1, 5) model.
In [172]:
perform_and_print_eacf <- function(da, ar.max, ma.max) {
    eacf_obj <- eacf(da, ar.max = ar.max, ma.max = ma.max)
    eacf_stats_tb = format(as.data.frame(eacf_obj$eacf), digits = 3)
    names(eacf_stats_tb) <- seq(from = 0, to = ma.max)
    display(eacf_stats_tb)
    display(eacf_obj$symbol)
    # pp67, asymptotic standard error of EACF
    display(2/sqrt(length(da)))
    c(eacf_obj, eacf_stats_tb)
}
In [173]:
unem_rate_eacf_res <- perform_and_print_eacf(diff_unem_rate, ar.max = 25, ma.max = 12)
A data.frame: 26 × 13
0123456789101112
<I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>>
1 0.1119 0.29465 0.21842 0.17198 1.94e-01 8.39e-02 0.05184 0.051878 0.02445-0.096187 0.026951-0.1594-0.047376
2-0.3443 0.15698-0.00881-0.03887 1.07e-01-3.19e-02-0.01686 0.022971 0.00774-0.077354-0.024625-0.1687 0.067674
3-0.4642 0.13045-0.01960-0.02493 1.13e-01-4.40e-02-0.06191 0.004640 0.01418-0.067424 0.002018-0.1414 0.040126
4-0.3650-0.35458 0.02198 0.08834 9.75e-02-5.66e-05 0.01257-0.004783 0.01615-0.066850-0.018419-0.1430-0.070736
5-0.4954-0.30539 0.04120 0.30176 9.52e-02-2.57e-03 0.01766 0.009470 0.04559-0.042195 0.047631-0.1245-0.039347
6 0.2517 0.44829 0.07516 0.16677 1.08e-01-8.10e-02-0.01934 0.000521 0.06292 0.022652-0.002533-0.1239-0.036224
7-0.3005 0.35966-0.13538-0.00859 3.02e-01-1.25e-01-0.01774 0.001050 0.05707 0.012168-0.005929-0.1340 0.021891
8-0.3330-0.00552-0.08958-0.02068 1.83e-01-1.10e-01 0.00728-0.032169 0.06110-0.009638 0.022461-0.0940 0.004250
9-0.3266 0.00149-0.33486-0.11095 1.03e-01 7.02e-02 0.24701 0.252584-0.02356-0.027958 0.056793-0.0928-0.000937
10-0.0757-0.16764-0.36087-0.22614 1.28e-01 2.94e-02 0.17613 0.275884-0.10130-0.033224 0.050234-0.0713-0.073394
11 0.2026-0.47690-0.28373 0.31492 1.15e-01-8.30e-02 0.01870 0.213078 0.19785 0.161355-0.028338-0.0807-0.073848
12 0.2286-0.47033-0.25334 0.19384 1.51e-01-1.63e-01 0.04894 0.212287 0.26369 0.145681-0.107458-0.0696 0.024317
13-0.1529 0.17761 0.20052 0.01742-2.30e-02 5.50e-02 0.11660-0.016648-0.10144 0.065135-0.049605-0.1695-0.024814
14 0.4946 0.10012 0.14184 0.13014 5.92e-05 1.60e-02 0.14602-0.014817-0.11259-0.022876-0.003969-0.1224 0.157803
15-0.3646 0.36006 0.08090 0.07057-2.23e-02 9.10e-03 0.01581 0.036602-0.11763 0.003957-0.005981-0.1541 0.218884
16-0.3713 0.26612-0.11562-0.00153 1.90e-02-1.23e-02 0.01686 0.018377-0.00921-0.015470 0.003599-0.1360 0.224311
17 0.3182 0.32908-0.11731-0.01477 8.46e-02-5.94e-02 0.01226-0.026056-0.00161-0.029129-0.007909-0.2376 0.189411
18-0.4914 0.15908-0.17106 0.04494 7.80e-02 1.66e-02 0.02160 0.014086-0.01989-0.000608-0.008598-0.3270-0.057128
19 0.2433 0.23853-0.15809 0.00850 3.95e-02 4.90e-02 0.00793 0.034214-0.02280-0.006053 0.000501-0.3291 0.077020
20-0.5005 0.42667-0.24634 0.00140 1.12e-01-1.23e-02 0.07434 0.057412-0.04946 0.004030 0.000927-0.2956 0.222547
21 0.4127 0.40983 0.09004 0.02077 5.28e-02 1.99e-02 0.06946 0.102419-0.03728-0.033084-0.063026-0.2526 0.120537
22-0.4559 0.18572-0.09269 0.00801 8.55e-02-4.82e-02 0.06485-0.051137 0.11544-0.129892-0.019547-0.2897-0.073481
23-0.4140-0.20397-0.03547 0.17919 1.34e-01-1.95e-01-0.07790 0.058271 0.21587-0.151503 0.016659-0.1912 0.124205
24-0.2072-0.26001-0.03256 0.22710 1.20e-01-1.92e-01-0.06062 0.140398 0.20160-0.156237-0.183888-0.2841 0.194636
25 0.0867 0.14692 0.37815 0.10849 8.49e-02 1.53e-01-0.04019-0.059815-0.11085-0.052183-0.286084-0.3222-0.139143
26-0.4406-0.00685 0.31666-0.11685-2.49e-02 1.73e-01-0.09294-0.008748-0.05978 0.004758-0.021080-0.3346 0.137259
A matrix: 26 × 13 of type chr
0123456789101112
0xxxxxxoooxoxo
1xxooxooooxoxo
2xxooxooooooxo
3xxoxxooooooxo
4xxoxxooooooxo
5xxxxxxoooooxo
6xxxoxxoooooxo
7xoxoxxoooooxo
8xoxxxoxxoooxo
9xxxxxoxxxoooo
10xxxxxxoxxxoxo
11xxxxxxoxxxxoo
12xxxoooxoxooxo
13xxxxooxoxooxx
14xxxoooooxooxx
15xxxooooooooxx
16xxxoxooooooxx
17xxxoxooooooxo
18xxxooooooooxx
19xxxoxooooooxx
20xxxooooxoooxx
21xxxoxoooxxoxo
22xxoxxxxoxxoxx
23xxoxxxoxxxxxx
24xxxxxxooxoxxx
25xoxxoxxooooxx
0.0738213470525562
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12
0  x x x x x x o o o x o  x  o 
1  x x o o x o o o o x o  x  o 
2  x x o o x o o o o o o  x  o 
3  x x o x x o o o o o o  x  o 
4  x x o x x o o o o o o  x  o 
5  x x x x x x o o o o o  x  o 
6  x x x o x x o o o o o  x  o 
7  x o x o x x o o o o o  x  o 
8  x o x x x o x x o o o  x  o 
9  x x x x x o x x x o o  o  o 
10 x x x x x x o x x x o  x  o 
11 x x x x x x o x x x x  o  o 
12 x x x o o o x o x o o  x  o 
13 x x x x o o x o x o o  x  x 
14 x x x o o o o o x o o  x  x 
15 x x x o o o o o o o o  x  x 
16 x x x o x o o o o o o  x  x 
17 x x x o x o o o o o o  x  o 
18 x x x o o o o o o o o  x  x 
19 x x x o x o o o o o o  x  x 
20 x x x o o o o x o o o  x  x 
21 x x x o x o o o x x o  x  o 
22 x x o x x x x o x x o  x  x 
23 x x o x x x o x x x x  x  x 
24 x x x x x x o o x o x  x  x 
25 x o x x o x x o o o o  x  x 

Fit an ARIMA(3,1,5) model

In [108]:
unem_rate_mod = stats::arima(unem_rate_ts, order = c(3,1,5), include.mean = F)
unem_rate_mod
Call:
stats::arima(x = unem_rate_ts, order = c(3, 1, 5), include.mean = F)

Coefficients:
        ar1     ar2      ar3      ma1      ma2     ma3     ma4     ma5
      0.671  0.7420  -0.7722  -0.6932  -0.5326  0.7862  -0.136  0.2068
s.e.  0.038  0.0289   0.0355   0.0505   0.0502  0.0493   0.049  0.0392

sigma^2 estimated as 0.03825:  log likelihood = 154.88,  aic = -291.77
In [110]:
unem_rate_mod_1 = arima(unem_rate_ts, order = c(3,1,5))
unem_rate_mod_1
Call:
arima(x = unem_rate_ts, order = c(3, 1, 5))

Coefficients:
        ar1     ar2      ar3      ma1      ma2     ma3     ma4     ma5
      0.671  0.7420  -0.7722  -0.6932  -0.5326  0.7862  -0.136  0.2068
s.e.  0.038  0.0289   0.0355   0.0505   0.0502  0.0493   0.049  0.0392

sigma^2 estimated as 0.03825:  log likelihood = 154.88,  aic = -293.77
In [111]:
diff_unem_rate_mod = arima(diff_unem_rate_ts, order = c(3,0,5), include.mean = F)
diff_unem_rate_mod
Call:
arima(x = diff_unem_rate_ts, order = c(3, 0, 5), include.mean = F)

Coefficients:
         ar1     ar2      ar3      ma1      ma2     ma3     ma4     ma5
      0.6711  0.7420  -0.7722  -0.6932  -0.5326  0.7862  -0.136  0.2068
s.e.  0.0379  0.0289   0.0355   0.0505   0.0502  0.0493   0.049  0.0392

sigma^2 estimated as 0.03825:  log likelihood = 154.88,  aic = -293.77
In [112]:
plot_time_fig(unem_rate_mod$residuals, main = "Residuals")
plot_acf(unem_rate_mod$residuals, main = "ACF of Residuals")
No description has been provided for this image
No description has been provided for this image

Box-Ljung test

  • pp33 very briefly talks about how to choose lag in the Box.test
  • From the following testing result, we cannot reject the null-hypothesis. So the model specified is adequate.
In [113]:
length(unem_rate); log(length(unem_rate))
735
6.59987049921284
In [114]:
Box.test(unem_rate_mod$residuals, lag = 7, type = 'Ljung')
	Box-Ljung test

data:  unem_rate_mod$residuals
X-squared = 5.5481, df = 7, p-value = 0.5934
In [115]:
Box.test(unem_rate_mod$residuals, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  unem_rate_mod$residuals
X-squared = 20.042, df = 12, p-value = 0.06629
  • Diagnose time-series ARIMA model
In [161]:
help(tsdiag)
In [162]:
par(bg = 'white')
tsdiag(unem_rate_mod, gof.lag=36)
No description has been provided for this image

Check business cycle

  • pp42
  • Compute the average length of business cycles
  • Doesn't matter if using the characteristic roots or the inverse of them.
  • Unit is "month". Average business cycle is about 14.2 months, which seems to be shoter than expected from the time plot above.
    • The business cycle is very sensitive to the order chosen in fitting ARIMA model. When I use ARIMA(2,0,5) or ARIMA(4,0,5), the business cycle doens't exist or too short.
In [69]:
find('arima')
  1. 'package:TSA'
  2. 'package:stats'
In [79]:
help("arima", package = "TSA")
In [80]:
help("arima", package = "stats")
In [116]:
sqrt(unem_rate_mod$sigma2)
0.195566624165521
In [117]:
unem_rate_mod$coef
ar1
0.671039253140843
ar2
0.741978426548034
ar3
-0.772184227253138
ma1
-0.693171209477596
ma2
-0.532634001385892
ma3
0.786216728216251
ma4
-0.135954862493427
ma5
0.206764447296786
In [118]:
ar_poly = c(1, -unem_rate_mod$coef[1:3]) # Characteristic equation for AR process
roots = polyroot(ar_poly)
roots
  1. 1.00437273858239+0.47656129169554i
  2. -1.04786282719477-0i
  3. 1.00437273858239-0.47656129169554i
In [119]:
for (i in 1:3) {
    print(paste("====", i))
    print(as.numeric(Mod(1-unem_rate_mod$coef[1]*roots[i]-unem_rate_mod$coef[2]*roots[i]^2-unem_rate_mod$coef[3]*roots[i]^3-unem_rate_mod$coef[4]*roots[i]^4)))
    print(as.numeric(Mod(roots[i])))
}
[1] "==== 1"
[1] 1.058741
[1] 1.111699
[1] "==== 2"
[1] 0.8357151
[1] 1.047863
[1] "==== 3"
[1] 1.058741
[1] 1.111699
  • Modulus of roots
In [120]:
roots; Mod(roots)
  1. 1.00437273858239+0.47656129169554i
  2. -1.04786282719477-0i
  3. 1.00437273858239-0.47656129169554i
  1. 1.11169926812516
  2. 1.04786282719477
  3. 1.11169926812516
In [121]:
1/roots; Mod(1/roots)
  1. 0.812681318944363-0.385606303531897i
  2. -9.54323384747883e-01+9e-16i
  3. 0.812681318944364+0.385606303531896i
  1. 0.899523844867207
  2. 0.954323384747883
  3. 0.899523844867207
  • Compute the average length of business cycles
  • Doesn't matter if using the characteristic roots or the inverse of them.
In [122]:
roots[1]; Re(roots[1]); Im(roots[1])
1.00437273858239+0.47656129169554i
1.00437273858239
0.476561291695542
In [125]:
k = 2*pi/acos(Re(roots[1])/Mod(roots[1]))
k
14.1823254232
In [126]:
k = 2*pi/acos(Re(1/roots[1])/Mod(1/roots[1]))
k
14.1823254232

Forecast

In [128]:
length(unem_rate_ts); unem_rate_ts
735
A Time Series: 62 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
1948 3.4 3.8 4.0 3.9 3.5 3.6 3.6 3.9 3.8 3.7 3.8 4.0
1949 4.3 4.7 5.0 5.3 6.1 6.2 6.7 6.8 6.6 7.9 6.4 6.6
1950 6.5 6.4 6.3 5.8 5.5 5.4 5.0 4.5 4.4 4.2 4.2 4.3
1951 3.7 3.4 3.4 3.1 3.0 3.2 3.1 3.1 3.3 3.5 3.5 3.1
1952 3.2 3.1 2.9 2.9 3.0 3.0 3.2 3.4 3.1 3.0 2.8 2.7
1953 2.9 2.6 2.6 2.7 2.5 2.5 2.6 2.7 2.9 3.1 3.5 4.5
1954 4.9 5.2 5.7 5.9 5.9 5.6 5.8 6.0 6.1 5.7 5.3 5.0
1955 4.9 4.7 4.6 4.7 4.3 4.2 4.0 4.2 4.1 4.3 4.2 4.2
1956 4.0 3.9 4.2 4.0 4.3 4.3 4.4 4.1 3.9 3.9 4.3 4.2
1957 4.2 3.9 3.7 3.9 4.1 4.3 4.2 4.1 4.4 4.5 5.1 5.2
1958 5.8 6.4 6.7 7.4 7.4 7.3 7.5 7.4 7.1 6.7 6.2 6.2
1959 6.0 5.9 5.6 5.2 5.1 5.0 5.1 5.2 5.5 5.7 5.8 5.3
1960 5.2 4.8 5.4 5.2 5.1 5.4 5.5 5.6 5.5 6.1 6.1 6.6
1961 6.6 6.9 6.9 7.0 7.1 6.9 7.0 6.6 6.7 6.5 6.1 6.0
1962 5.8 5.5 5.6 5.6 5.5 5.5 5.4 5.7 5.6 5.4 5.7 5.5
1963 5.7 5.9 5.7 5.7 5.9 5.6 5.6 5.4 5.5 5.5 5.7 5.5
1964 5.6 5.4 5.4 5.3 5.1 5.2 4.9 5.0 5.1 5.1 4.8 5.0
1965 4.9 5.1 4.7 4.8 4.6 4.6 4.4 4.4 4.3 4.2 4.1 4.0
1966 4.0 3.8 3.8 3.8 3.9 3.8 3.8 3.8 3.7 3.7 3.6 3.8
1967 3.9 3.8 3.8 3.8 3.8 3.9 3.8 3.8 3.8 4.0 3.9 3.8
1968 3.7 3.8 3.7 3.5 3.5 3.7 3.7 3.5 3.4 3.4 3.4 3.4
1969 3.4 3.4 3.4 3.4 3.4 3.5 3.5 3.5 3.7 3.7 3.5 3.5
1970 3.9 4.2 4.4 4.6 4.8 4.9 5.0 5.1 5.4 5.5 5.9 6.1
1971 5.9 5.9 6.0 5.9 5.9 5.9 6.0 6.1 6.0 5.8 6.0 6.0
1972 5.8 5.7 5.8 5.7 5.7 5.7 5.6 5.6 5.5 5.6 5.3 5.2
1973 4.9 5.0 4.9 5.0 4.9 4.9 4.8 4.8 4.8 4.6 4.8 4.9
1974 5.1 5.2 5.1 5.1 5.1 5.4 5.5 5.5 5.9 6.0 6.6 7.2
1975 8.1 8.1 8.6 8.8 9.0 8.8 8.6 8.4 8.4 8.4 8.3 8.2
1976 7.9 7.7 7.6 7.7 7.4 7.6 7.8 7.8 7.6 7.7 7.8 7.8
1977 7.5 7.6 7.4 7.2 7.0 7.2 6.9 7.0 6.8 6.8 6.8 6.4
1978 6.4 6.3 6.3 6.1 6.0 5.9 6.2 5.9 6.0 5.8 5.9 6.0
1979 5.9 5.9 5.8 5.8 5.6 5.7 5.7 6.0 5.9 6.0 5.9 6.0
1980 6.3 6.3 6.3 6.9 7.5 7.6 7.8 7.7 7.5 7.5 7.5 7.2
1981 7.5 7.4 7.4 7.2 7.5 7.5 7.2 7.4 7.6 7.9 8.3 8.5
1982 8.6 8.9 9.0 9.3 9.4 9.6 9.8 9.810.110.410.810.8
198310.410.410.310.210.110.1 9.4 9.5 9.2 8.8 8.5 8.3
1984 8.0 7.8 7.8 7.7 7.4 7.2 7.5 7.5 7.3 7.4 7.2 7.3
1985 7.3 7.2 7.2 7.3 7.2 7.4 7.4 7.1 7.1 7.1 7.0 7.0
1986 6.7 7.2 7.2 7.1 7.2 7.2 7.0 6.9 7.0 7.0 6.9 6.6
1987 6.6 6.6 6.6 6.3 6.3 6.2 6.1 6.0 5.9 6.0 5.8 5.7
1988 5.7 5.7 5.7 5.4 5.6 5.4 5.4 5.6 5.4 5.4 5.3 5.3
1989 5.4 5.2 5.0 5.2 5.2 5.3 5.2 5.2 5.3 5.3 5.4 5.4
1990 5.4 5.3 5.2 5.4 5.4 5.2 5.5 5.7 5.9 5.9 6.2 6.3
1991 6.4 6.6 6.8 6.7 6.9 6.9 6.8 6.9 6.9 7.0 7.0 7.3
1992 7.3 7.4 7.4 7.4 7.6 7.8 7.7 7.6 7.6 7.3 7.4 7.4
1993 7.3 7.1 7.0 7.1 7.1 7.0 6.9 6.8 6.7 6.8 6.6 6.5
1994 6.6 6.6 6.5 6.4 6.1 6.1 6.1 6.0 5.9 5.8 5.6 5.5
1995 5.6 5.4 5.4 5.8 5.6 5.6 5.7 5.7 5.6 5.5 5.6 5.6
1996 5.6 5.5 5.5 5.6 5.6 5.3 5.5 5.1 5.2 5.2 5.4 5.4
1997 5.3 5.2 5.2 5.1 4.9 5.0 4.9 4.8 4.9 4.7 4.6 4.7
1998 4.6 4.6 4.7 4.3 4.4 4.5 4.5 4.5 4.6 4.5 4.4 4.4
1999 4.3 4.4 4.2 4.3 4.2 4.3 4.3 4.2 4.2 4.1 4.1 4.0
2000 4.0 4.1 4.0 3.8 4.0 4.0 4.0 4.1 3.9 3.9 3.9 3.9
2001 4.2 4.2 4.3 4.4 4.3 4.5 4.6 4.9 5.0 5.3 5.5 5.7
2002 5.7 5.7 5.7 5.9 5.8 5.8 5.8 5.7 5.7 5.7 5.9 6.0
2003 5.8 5.9 5.9 6.0 6.1 6.3 6.2 6.1 6.1 6.0 5.8 5.7
2004 5.7 5.6 5.8 5.6 5.6 5.6 5.5 5.4 5.4 5.5 5.4 5.4
2005 5.2 5.4 5.2 5.2 5.1 5.1 5.0 4.9 5.0 5.0 5.0 4.8
2006 4.7 4.8 4.7 4.7 4.7 4.6 4.7 4.7 4.5 4.4 4.5 4.4
2007 4.6 4.5 4.4 4.5 4.4 4.6 4.6 4.6 4.7 4.7 4.7 5.0
2008 5.0 4.8 5.1 5.0 5.4 5.5 5.8 6.1 6.2 6.6 6.9 7.4
2009 7.7 8.2 8.6
In [140]:
npts = 8
eotr = length(unem_rate_ts)
h = 4
freq = 12
order = c(3,1,5)
fixed = NULL
seasonal = list(order = c(0L, 0L, 0L), period = NA)
unem_rate_fc_res = plot_arima_forecast_fig(
    da_ts=unem_rate_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', 
    include.mean=F, transform.pars=NULL,
    main="Forecasts from ARIMA(3,1,5) for unemp_rate", 
    xlab="Year", ylab="Unemp-Rate", ylim=c(5, 11)
)
[1] 734
2008.5 ; 2009.583
No description has been provided for this image

Alternative solutions (only use AR model)

In [147]:
par(bg = 'white')
# decay slowly
acf(unem_rate)
par(mfrow = c(2, 1))
acf(diff(unem_rate))
pacf(diff(unem_rate))
acf(diff(unem_rate, 12))
pacf(diff(unem_rate, 12))
cat("order =", ar(unem_rate, method = 'mle')$order, "\n")
m1 = arima(unem_rate_ts, order = c(11, 0, 0))
m1
# Follow the pp51 to modify and refine the model based on previous fit
m1 = arima(unem_rate_ts, order = c(11, 0, 0), fixed = c(NA, NA, 0, 0, 0, NA, 0, 0, 0, NA, NA, NA))
m1
# FIXME
m2 = arima(unem_rate_ts, order = c(2, 1, 1), seasonal = list(order = c(1, 0, 1), period = 12))
m2
tsdiag(m1, gof = 36)
tsdiag(m2, gof = 36)
predict(m1, 4)
predict(m2, 4)
No description has been provided for this image
No description has been provided for this image
order = 11 
Call:
arima(x = unem_rate_ts, order = c(11, 0, 0))

Coefficients:
         ar1     ar2      ar3      ar4     ar5      ar6      ar7     ar8
      0.9886  0.2375  -0.0741  -0.0630  0.0301  -0.1283  -0.0426  0.0539
s.e.  0.0367  0.0516   0.0525   0.0525  0.0526   0.0524   0.0526  0.0527
          ar9     ar10    ar11  intercept
      -0.0146  -0.1293  0.1259     5.6554
s.e.   0.0526   0.0518  0.0371     0.4422

sigma^2 estimated as 0.03867:  log likelihood = 150.03,  aic = -276.07
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = unem_rate_ts, order = c(11, 0, 0), fixed = c(NA, NA, 0, 0, 0, NA, 
    0, 0, 0, NA, NA, NA))

Coefficients:
         ar1     ar2  ar3  ar4  ar5      ar6  ar7  ar8  ar9     ar10    ar11
      0.9805  0.1695    0    0    0  -0.1728    0    0    0  -0.1216  0.1282
s.e.  0.0361  0.0432    0    0    0   0.0270    0    0    0   0.0433  0.0364
      intercept
         5.6605
s.e.     0.4337

sigma^2 estimated as 0.03904:  log likelihood = 146.55,  aic = -281.11
Call:
arima(x = unem_rate_ts, order = c(2, 1, 1), seasonal = list(order = c(1, 0, 
    1), period = 12))

Coefficients:
         ar1     ar2      ma1    sar1     sma1
      0.5801  0.2444  -0.5813  0.5609  -0.8185
s.e.  0.0637  0.0398   0.0588  0.0734   0.0538

sigma^2 estimated as 0.03725:  log likelihood = 164.28,  aic = -318.56
No description has been provided for this image
No description has been provided for this image
$pred
A Time Series: 1 × 4
AprMayJunJul
20098.7967658.9819969.1125559.246428
$se
A Time Series: 1 × 4
AprMayJunJul
20090.19758560.27671340.35565800.4358310
$pred
A Time Series: 1 × 4
AprMayJunJul
20098.8418449.0171549.1269179.231977
$se
A Time Series: 1 × 4
AprMayJunJul
20090.19300980.27280390.36323950.4508560
No description has been provided for this image
In [149]:
ar_poly_1 = c(1, -m1$coef[1:11]) # Characteristic equation for AR process
roots_1 = polyroot(ar_poly_1)
roots_1
  1. 0.72130962198895+1.08701346704084i
  2. -1.14274870161967+0.39895294351557i
  3. -0.772376200376323-0.973913718979966i
  4. 1.135449573239+0.20175780984644i
  5. 0.01445275970928+1.21215403312275i
  6. -0.772376200376323+0.973913718979966i
  7. 0.01445275970928-1.21215403312275i
  8. 1.03577871789383+0i
  9. 1.135449573239-0.20175780984644i
  10. -1.14274870161967-0.39895294351557i
  11. 0.72130962198895-1.08701346704084i
In [155]:
Im(roots_1[8])
3.31804353375726e-15
In [156]:
for (i in 1:length(roots_1)) {
    if (abs(Im(roots_1[i]))>1e-8) {
        print(i)
        print(roots_1[i])
        print(2*pi/acos(Re(roots_1[i])/Mod(roots_1[i])))
    }
}
[1] 1
[1] 0.72131+1.087013i
[1] 6.379253
[1] 2
[1] -1.142749+0.398953i
[1] 2.239432
[1] 3
[1] -0.7723762-0.9739137i
[1] 2.803374
[1] 4
[1] 1.13545+0.201758i
[1] 35.72949
[1] 5
[1] 0.014453+1.212154i
[1] 4.030593
[1] 6
[1] -0.7723762+0.9739137i
[1] 2.803374
[1] 7
[1] 0.014453-1.212154i
[1] 4.030593
[1] 9
[1] 1.13545-0.201758i
[1] 35.72949
[1] 10
[1] -1.142749-0.398953i
[1] 2.239432
[1] 11
[1] 0.72131-1.087013i
[1] 6.379253
In [151]:
ar_poly_2 = c(1, -m2$coef[1:2]) # Characteristic equation for AR process
roots_2 = polyroot(ar_poly_2)
roots_2
  1. 1.158427203051+0i
  2. -3.53236324765541-0i

2-4

In [3]:
da = read.table("../AFTS_data/Ch02/m-deciles08.txt", header = T)
dim(da); da[1:5,]
  1. 468
  2. 5
A data.frame: 5 × 5
dateCAP1RETCAP2RETCAP9RETCAP10RET
<int><dbl><dbl><dbl><dbl>
119700130 0.054383-0.004338-0.073082-0.076874
219700227 0.020264 0.020155 0.064185 0.059512
319700331-0.031790-0.028090-0.004034-0.001327
419700430-0.184775-0.193004-0.115825-0.091112
519700529-0.088189-0.085342-0.085565-0.053193
In [4]:
dec2_ts = ts(da$CAP2RET, frequency = 12, start = c(1970, 1))
dec10_ts = ts(da$CAP10RET, frequency = 12, start = c(1970, 1))
In [5]:
plot_time_fig(dec2_ts, main = "CRSP Decile Index", xlab = "Year")
lines(dec10_ts, type = , col = 'red', lty = 2)
legend(
    "topright", 
    legend = c("Decile 2", "Decile 10"), 
    col = c("black", "red"), 
    lty = c(1, 2),
    pch = c(8, NA)
)
No description has been provided for this image

(a) Test autocorrelations

  • Reject the null-hypothesis (i.e., the first 12 lags have no serial correlation) for Decile 2. So for Decile 2, the first 12 lags have serial correlation.
  • Cannot reject the null-hypothesis for Decile 10. So for Decile 10, the first 12 lags have no significant serial correlation.
In [6]:
par(mfrow = c(2, 1), bg = 'white')
acf(da$CAP2RET)
pacf(da$CAP2RET)
Box.test(da$CAP2RET, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  da$CAP2RET
X-squared = 55.736, df = 12, p-value = 1.335e-07
No description has been provided for this image
In [7]:
par(mfrow = c(2, 1), bg = 'white')
acf(da$CAP10RET)
pacf(da$CAP10RET)
Box.test(da$CAP10RET, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  da$CAP10RET
X-squared = 10.687, df = 12, p-value = 0.5559
No description has been provided for this image

(b) Fit a model and perform model checking

In [8]:
dec2_eacf_res <- perform_and_print_eacf(da$CAP2RET, ar.max = 25, ma.max = 12)
A data.frame: 26 × 13
0123456789101112
<I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>>
1 0.20369-0.03114-0.0617-0.03854-0.04448-0.01246-1.80e-02-0.08682-0.05043 0.044819 0.03371 0.23095 0.02117
2 0.32892-0.10631-0.0375 0.01921-0.03447 0.02160 3.98e-03-0.06861-0.05926 0.066053-0.00460 0.23682 0.06396
3-0.42524 0.01475-0.1102-0.01244-0.01521 0.00803 7.46e-03-0.04457-0.04307 0.020265 0.01009 0.23522-0.05571
4-0.40167-0.05250-0.0518 0.06362-0.00120 0.01725 2.15e-03-0.02434-0.00964 0.024374-0.02813 0.23364-0.13191
5-0.40022-0.43254-0.2838-0.00515-0.03041 0.01686 2.73e-03-0.00353-0.02382-0.007906-0.00472 0.23347-0.06624
6-0.06468-0.43194-0.1068 0.12936-0.06677-0.04765-4.82e-03 0.00172-0.02487 0.000663-0.00560 0.23064-0.03428
7-0.10539-0.02918-0.3781-0.13146-0.13792 0.01120-1.39e-02-0.00394-0.02474-0.006610-0.00960 0.21633-0.12140
8-0.25901 0.04140-0.3007 0.12273-0.00148-0.11108-4.85e-05 0.03159 0.00131-0.003169-0.02814 0.17372-0.12466
9-0.24123 0.39312-0.2237 0.09617-0.00180-0.10194 6.02e-05 0.02881 0.01067-0.003566-0.03400 0.15280-0.10936
10 0.39313 0.35825 0.0508 0.25484 0.06186-0.05278 1.18e-01 0.03313 0.00458-0.005226-0.04365 0.08590-0.12888
11 0.02324-0.16438 0.0591 0.29471 0.14351-0.01095 1.10e-01 0.03473-0.04732 0.039074-0.02754 0.01843-0.12941
12 0.00277-0.19923 0.1018 0.30826 0.02356-0.05258 1.38e-01 0.12985 0.08072 0.153820-0.42198-0.00799-0.06233
13 0.29698 0.23296 0.0675 0.05170-0.05332 0.08408 2.25e-01 0.12233 0.20471 0.187092 0.07290 0.02734-0.08006
14-0.41434 0.15559-0.1088 0.03756 0.00361 0.02920 2.07e-01-0.10366 0.03855 0.207229 0.03509-0.18728 0.01939
15 0.27098-0.00392 0.0464 0.00807-0.01820 0.03644 1.96e-01-0.10690 0.03093 0.170090 0.04863-0.03330 0.01009
16 0.32328-0.00167 0.0780-0.00315-0.08522 0.00227 1.16e-01 0.06426 0.10021 0.061636 0.05065-0.04054-0.02312
17 0.19872-0.24388 0.0253-0.03706-0.08837 0.00731 5.98e-02-0.03659 0.08578-0.027778 0.11103-0.01494-0.01742
18 0.45276 0.05983 0.0949-0.08135-0.26125 0.03599 1.02e-01 0.01364 0.05305 0.031292 0.11314-0.07872-0.00302
19-0.39023 0.35687-0.2743 0.26657-0.15358 0.05311 9.46e-02-0.10945-0.01772 0.034096 0.13725 0.02100 0.00207
20 0.13505-0.37359-0.4908 0.21703 0.27756 0.18885 2.24e-03-0.06054-0.04551 0.019328 0.07323-0.08491 0.01630
21 0.21470-0.18877-0.4752 0.14720 0.39244 0.15284-5.28e-03-0.12773 0.03177-0.008420 0.07184-0.00911 0.00216
22-0.46086 0.30565-0.1895 0.15917 0.36093-0.05730 1.77e-01-0.08461 0.02247-0.002969 0.04152-0.04145-0.02601
23 0.37876 0.48961 0.3502 0.07946 0.13617-0.10137 2.15e-01-0.05946 0.00999 0.043435 0.06366-0.01055-0.01546
24-0.37059 0.39806-0.1081 0.19561 0.26711 0.10149 2.48e-01 0.00868 0.10901 0.107988 0.06787-0.04747 0.07475
25-0.22999-0.42493 0.0483 0.38694-0.16429-0.24948 2.03e-01 0.28052 0.08915-0.240086-0.00779 0.21276-0.04742
26-0.16376-0.37835 0.0813 0.30938-0.14831-0.34927 1.43e-01 0.26769 0.06260-0.239522-0.01508 0.10886-0.08528
A matrix: 26 × 13 of type chr
0123456789101112
0xooooooooooxo
1xxoooooooooxo
2xoxooooooooxo
3xooooooooooxx
4xxxooooooooxo
5oxxxoooooooxo
6xoxxxooooooxx
7xoxxoxoooooxx
8xxxxoxoooooxx
9xxoxooxooooox
10oxoxxoxooooox
11oxxxooxxoxxoo
12xxooooxxxxooo
13xxxoooxxoxoxo
14xoooooxxoxooo
15xoooooxoxoooo
16xxooooooooxoo
17xoxoxoxoooxoo
18xxxxxooxooxoo
19xxxxxxooooooo
20xxxxxxoxooooo
21xxxxxoxoooooo
22xxxoxxxoooooo
23xxxxxxxoxxooo
24xxoxxxxxoxoxo
25xxoxxxxxoxoxo
0.0924500327042048
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12
0  x o o o o o o o o o o  x  o 
1  x x o o o o o o o o o  x  o 
2  x o x o o o o o o o o  x  o 
3  x o o o o o o o o o o  x  x 
4  x x x o o o o o o o o  x  o 
5  o x x x o o o o o o o  x  o 
6  x o x x x o o o o o o  x  x 
7  x o x x o x o o o o o  x  x 
8  x x x x o x o o o o o  x  x 
9  x x o x o o x o o o o  o  x 
10 o x o x x o x o o o o  o  x 
11 o x x x o o x x o x x  o  o 
12 x x o o o o x x x x o  o  o 
13 x x x o o o x x o x o  x  o 
14 x o o o o o x x o x o  o  o 
15 x o o o o o x o x o o  o  o 
16 x x o o o o o o o o x  o  o 
17 x o x o x o x o o o x  o  o 
18 x x x x x o o x o o x  o  o 
19 x x x x x x o o o o o  o  o 
20 x x x x x x o x o o o  o  o 
21 x x x x x o x o o o o  o  o 
22 x x x o x x x o o o o  o  o 
23 x x x x x x x o x x o  o  o 
24 x x o x x x x x o x o  x  o 
25 x x o x x x x x o x o  x  o 
  • Fit ARIMA(1,0,1) model, b/c only this model has all coefs significant.
In [10]:
dec2_mod_1 <- arima(da$CAP2RET, order = c(1, 0, 1), include.mean = T)
dec2_mod_1; Box.test(dec2_mod_1$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_1, gof.lag = 36)
Call:
arima(x = da$CAP2RET, order = c(1, 0, 1), include.mean = T)

Coefficients:
          ar1     ma1  intercept
      -0.0260  0.2433     0.0105
s.e.   0.1773  0.1702     0.0036

sigma^2 estimated as 0.004077:  log likelihood = 623.48,  aic = -1240.96
	Box-Ljung test

data:  dec2_mod_1$residuals
X-squared = 36.405, df = 12, p-value = 0.0002788
No description has been provided for this image
In [11]:
# Follow the pp51 to modify and refine the model based on previous fit
dec2_mod_2 <- arima(da$CAP2RET, order = c(0, 0, 1), include.mean = T)
dec2_mod_2; Box.test(dec2_mod_2$residuals, lag = 10, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_2, gof.lag = 36)
Call:
arima(x = da$CAP2RET, order = c(0, 0, 1), include.mean = T)

Coefficients:
         ma1  intercept
      0.2191     0.0105
s.e.  0.0449     0.0036

sigma^2 estimated as 0.004077:  log likelihood = 623.47,  aic = -1242.93
	Box-Ljung test

data:  dec2_mod_2$residuals
X-squared = 8.1094, df = 10, p-value = 0.6182
No description has been provided for this image
In [12]:
dec2_mod_3 <- arima(da$CAP2RET, order = c(4, 0, 3), include.mean = T)
dec2_mod_3; Box.test(dec2_mod_3$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_3, gof.lag = 36)
Call:
arima(x = da$CAP2RET, order = c(4, 0, 3), include.mean = T)

Coefficients:
         ar1     ar2      ar3     ar4      ma1      ma2     ma3  intercept
      0.7561  0.7891  -0.7599  0.2034  -0.5337  -0.9888  0.5226     0.0117
s.e.  0.6190  0.1121   0.6092  0.1197   0.6306   0.0121  0.6263     0.0013

sigma^2 estimated as 0.003948:  log likelihood = 629.86,  aic = -1243.71
	Box-Ljung test

data:  dec2_mod_3$residuals
X-squared = 30.035, df = 12, p-value = 0.002759
No description has been provided for this image
In [13]:
dec2_mod_3 <- arima(
    da$CAP2RET, 
    order = c(4, 0, 3), 
    # Follow the pp51 to modify and refine the model based on previous fit
    fixed = c(0, NA, 0, NA, 0, NA, 0, NA),
    include.mean = T
)
dec2_mod_3; Box.test(dec2_mod_3$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_3, gof.lag = 36)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = da$CAP2RET, order = c(4, 0, 3), include.mean = T, fixed = c(0, NA, 
    0, NA, 0, NA, 0, NA))

Coefficients:
      ar1     ar2  ar3      ar4  ma1      ma2  ma3  intercept
        0  0.3407    0  -0.0362    0  -0.3752    0     0.0106
s.e.    0  0.3787    0   0.0507    0   0.3767    0     0.0027

sigma^2 estimated as 0.004255:  log likelihood = 613.47,  aic = -1218.94
	Box-Ljung test

data:  dec2_mod_3$residuals
X-squared = 51.445, df = 12, p-value = 7.773e-07
No description has been provided for this image

(c) Forecast using ARMA model

In [22]:
npts = 12
eotr = length(dec2_ts)-npts
h = npts
freq = 12
order = c(0,0,1)
fixed = NULL
seasonal = NULL
dec2_fc_res = plot_arima_forecast_fig(
    da_ts=dec2_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', 
    include.mean=T, transform.pars=NULL,
    main="Forecasts from ARIMA(0,0,1) for Decile 2 index", 
    xlab="Year", ylab="Unemp-Rate"# , ylim=c(5, 11)
)
dec2_fc_tb = comb_forecast_res(dec2_fc_res, dec2_ts, eotr, freq)
dec2_fc_tb
[1] 456
2006.917 ; 2009
Forecast method: ARIMA(0,0,1) with non-zero mean

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
         ma1  intercept
      0.2043     0.0119
s.e.  0.0464     0.0036

sigma^2 estimated as 0.004074:  log likelihood = 607.66,  aic = -1211.33

Error measures:
                       ME       RMSE       MAE      MPE     MAPE      MASE
Training set 1.545115e-06 0.06382686 0.0433542 100.5658 186.3821 0.7337148
                     ACF1
Training set -0.004704323

Forecasts:
         Point Forecast       Lo 80      Hi 80      Lo 95     Hi 95
Jan 2008    0.009250982 -0.07254643 0.09104840 -0.1158474 0.1343493
Feb 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Mar 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Apr 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
May 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Jun 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Jul 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Aug 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Sep 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Oct 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Nov 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
Dec 2008    0.011867122 -0.07161973 0.09535397 -0.1158150 0.1395492
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.0092509820.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.011867122
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.063826860.065145130.065145130.065145130.065145130.065145130.065145130.065145130.065145130.065145130.065145130.06514513
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
2008-0.008644-0.017025-0.024443 0.015660 0.034079-0.037152-0.048765-0.012646-0.132327-0.129844-0.114484-0.022210
A Time Series: 12 × 3
ForecastStd. ErrorActual
Jan 20080.0092509820.06382686-0.008644
Feb 20080.0118671220.06514513-0.017025
Mar 20080.0118671220.06514513-0.024443
Apr 20080.0118671220.06514513 0.015660
May 20080.0118671220.06514513 0.034079
Jun 20080.0118671220.06514513-0.037152
Jul 20080.0118671220.06514513-0.048765
Aug 20080.0118671220.06514513-0.012646
Sep 20080.0118671220.06514513-0.132327
Oct 20080.0118671220.06514513-0.129844
Nov 20080.0118671220.06514513-0.114484
Dec 20080.0118671220.06514513-0.022210
No description has been provided for this image
In [17]:
class(dec2_fc_res); attributes(dec2_fc_res); dec2_fc_res$model
'forecast'
$names
  1. 'method'
  2. 'model'
  3. 'level'
  4. 'mean'
  5. 'lower'
  6. 'upper'
  7. 'x'
  8. 'series'
  9. 'fitted'
  10. 'residuals'
$class
'forecast'
Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
         ma1  intercept
      0.2043     0.0119
s.e.  0.0464     0.0036

sigma^2 estimated as 0.004074:  log likelihood = 607.66,  aic = -1211.33
In [20]:
help("predict")
In [21]:
predict(dec2_fc_res$model, 12)
$pred
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.0092509820.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.0118671220.011867122
$se
A Time Series: 1 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20080.063826860.065145130.065145130.065145130.065145130.065145130.065145130.065145130.065145130.065145130.065145130.06514513
In [23]:
dec2_fc_tb
A Time Series: 12 × 3
ForecastStd. ErrorActual
Jan 20080.0092509820.06382686-0.008644
Feb 20080.0118671220.06514513-0.017025
Mar 20080.0118671220.06514513-0.024443
Apr 20080.0118671220.06514513 0.015660
May 20080.0118671220.06514513 0.034079
Jun 20080.0118671220.06514513-0.037152
Jul 20080.0118671220.06514513-0.048765
Aug 20080.0118671220.06514513-0.012646
Sep 20080.0118671220.06514513-0.132327
Oct 20080.0118671220.06514513-0.129844
Nov 20080.0118671220.06514513-0.114484
Dec 20080.0118671220.06514513-0.022210

2-5

  • There is some long-range dependency.
In [25]:
da = read.table("../AFTS_sol/data/d-ibm3dx7008.txt", header = T)
ibm = da$rtn
plot_acf(abs(ibm), lag.max = 100)
No description has been provided for this image

2-6

  • This is a multiplicative seasonal model.
  • The seasonal prediction is very good.
  • The data set in this problem bears a lot similarity with the JnJ quarterly earning data set in Sec. 2.8. After taking one regular difference, the seasonality becomes very obvious. This means that after taking one regular difference, we need to model the AR or MA pattern of the resulting time series.
  • The logic of determining the orders in arima model from the various pacf and acf plots can be as following.
    • The plots for pow show strong serial correlation. So we need to take diff.
    • The plots for diff(pow) show strong seasonality in the acf. So we need to take the seasonal diff.
    • The plots for diff(pow, 12) show some autocorrelation in the acf. So we need to take both diff and seasonal diff.
    • The plots for diff(diff(pow), 12) show some autocorrelation in the acf. So we need to apply MA model on diff(diff(pow), 12). Thus this is a multiplicative seasonal model.
In [100]:
da = read.table("../AFTS_sol/data/power6.txt", header = F)
dim(da); da[1:5,]
  1. 264
  2. 1
  1. 9.2691
  2. 9.2465
  3. 9.2639
  4. 9.2784
  5. 9.3147
In [101]:
par(bg = 'white')
pow = da[, 1]
pow_ts = ts(pow, frequency = 12, start = c(2000, 1))
# From the time plot, should immediately think that we can apply multiplicative seasonal model
plot(pow_ts, type = 'o', xlab = 'Year', ylab = 'Power', main = "Power Consumption")

plot(diff(pow_ts), type = 'o', xlab = 'Year', ylab = 'diff(Power)', main = "Power Consumption")

plot(diff(pow_ts, 12), type = 'o', xlab = 'Year', ylab = 'diff(Power, 12)', main = "Power Consumption")

plot(diff(diff(pow_ts), 12), type = 'o', xlab = 'Year', ylab = 'diff(diff(Power), 12)', main = "Power Consumption")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [102]:
par(mfrow = c(2, 1), bg = 'white')
pacf(pow)
acf(pow)
pacf(diff(pow))
acf(diff(pow))
pacf(diff(pow, 12))
acf(diff(pow, 12))
pacf(diff(diff(pow), 12))
acf(diff(diff(pow), 12))
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [103]:
# The eacf is not helpful in determining the order of ARIMA model
pow_eacf_res <- perform_and_print_eacf(pow, ar.max = 25, ma.max = 12)
A data.frame: 26 × 13
0123456789101112
<I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>>
1 0.96571 0.9080 0.8377 0.7655 0.71498 0.69168 0.69440 0.72362 0.773267 0.8211 0.8573 0.86632 0.8387
2 0.48462 0.2868 0.0751-0.3656-0.49156-0.47986-0.42300-0.27176 0.055811 0.2656 0.5053 0.75737 0.5252
3-0.12526 0.2006 0.0500-0.2874-0.08769-0.01527-0.13087-0.29585 0.084182 0.0655 0.0642 0.48727 0.2065
4 0.40583 0.0580 0.1825-0.2774-0.05706-0.00408-0.00607-0.22609 0.164512-0.0636-0.0280 0.40705-0.0743
5-0.23966 0.1946 0.5000-0.2768-0.04448 0.05621-0.00483-0.21901 0.179330-0.0263-0.0688 0.41287-0.0364
6-0.40472-0.1180-0.3035-0.1740 0.09863-0.05222 0.10280-0.23891-0.089101-0.0882 0.0255 0.25582 0.0064
7-0.49066 0.0998-0.0943-0.1726 0.32596 0.16719 0.15146-0.22323 0.052400-0.1301 0.0103 0.24534 0.0243
8-0.43350-0.1371-0.0570-0.2257 0.39052-0.02031 0.10420-0.21595-0.124124-0.0587-0.0169 0.15903-0.0706
9-0.46920 0.0606-0.1957-0.2735 0.39116-0.01460 0.29364-0.19086 0.088586-0.0716-0.0433 0.22520-0.0638
10-0.40164-0.1888 0.1925 0.4526-0.14736-0.19761 0.16804 0.23199-0.063995-0.1397 0.0650-0.02067 0.0570
11-0.50251-0.3665 0.0424 0.5069-0.18772-0.32381 0.12976 0.28136-0.066008-0.2720 0.0595-0.00967-0.0380
12-0.01093 0.1967-0.1065 0.3678 0.14461 0.12197 0.07902 0.00539 0.137283-0.0978 0.0701 0.06164-0.0754
13-0.17885 0.2982 0.2797 0.3059-0.11268 0.03445 0.05553-0.02337 0.112943-0.1702 0.1513 0.03375-0.0246
14-0.49983 0.2406-0.2440 0.1704-0.10629 0.04095 0.07715-0.02402 0.032608-0.0461 0.0347-0.22437 0.1587
15-0.22348-0.3658-0.2354-0.0464 0.04177 0.08093 0.10552 0.09814-0.003118-0.0311-0.0309-0.38593 0.0145
16-0.40200-0.2337-0.1751-0.1385 0.02525-0.01050 0.11290 0.22483-0.066749-0.0521-0.0354-0.39753 0.1508
17-0.47804 0.1087-0.0335 0.0536 0.01371-0.00288 0.12952 0.20014-0.026307 0.0209-0.0252-0.32491 0.0708
18-0.43757-0.0678 0.0226 0.0106 0.00705-0.02970 0.16429 0.02475-0.031768 0.0287-0.0378-0.37439 0.1403
19-0.44875-0.0350 0.0283 0.0140 0.04629-0.16741-0.00992 0.03302 0.051821-0.0152-0.2024-0.40911-0.1789
20-0.10704 0.0546-0.0268 0.0597 0.26460-0.17678 0.01269 0.01676 0.053342-0.0303-0.1867-0.37761 0.1276
21 0.00603 0.1789 0.1906 0.3505 0.09617-0.05645 0.08630-0.05548 0.005178-0.0115-0.1248-0.33770 0.0559
22-0.34425-0.4013 0.1859 0.4218 0.14376-0.04038 0.21477-0.06536 0.000518-0.0341-0.0703-0.34027-0.0386
23-0.40302-0.4847 0.0806 0.4274 0.04296-0.15622 0.10458-0.07416-0.106884-0.0618-0.0455-0.34531-0.0273
24-0.24496 0.3079-0.1119 0.4502 0.14402 0.19132 0.11027-0.07358-0.103901-0.0729 0.0061-0.34701-0.1540
25 0.09975 0.4002 0.1389 0.4044-0.13304 0.10269 0.00357-0.08368 0.014541-0.2153-0.1404-0.35439 0.1572
26-0.44069-0.0591-0.1255 0.3419-0.20926 0.10856 0.00282-0.03992 0.104176-0.1713 0.1495-0.32197-0.1242
A matrix: 26 × 13 of type chr
0123456789101112
0xxxxxxxxxxxxx
1xxoxxxxxoxxxx
2xxoxooxxoooxx
3xoxxoooxxooxo
4xxxxoooxxooxo
5xoxxoooxoooxo
6xooxxxxxoxoxo
7xxoxxooxoooxo
8xoxxxoxxoooxo
9xxxxxxxxoxooo
10xxoxxxxxoxooo
11oxoxxoooxoooo
12xxxxoooooxxoo
13xxxxoooooooxx
14xxxooooooooxo
15xxxxoooxoooxx
16xoooooxxoooxo
17xoooooxooooxx
18xooooxooooxxx
19ooooxxooooxxo
20oxxxoooooooxo
21xxxxxoxooooxo
22xxoxoxoooooxo
23xxoxxxoooooxx
24oxxxxooooxxxx
25xooxxooooxxxo
0.123091490979333
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12
0  x x x x x x x x x x x  x  x 
1  x x o x x x x x o x x  x  x 
2  x x o x o o x x o o o  x  x 
3  x o x x o o o x x o o  x  o 
4  x x x x o o o x x o o  x  o 
5  x o x x o o o x o o o  x  o 
6  x o o x x x x x o x o  x  o 
7  x x o x x o o x o o o  x  o 
8  x o x x x o x x o o o  x  o 
9  x x x x x x x x o x o  o  o 
10 x x o x x x x x o x o  o  o 
11 o x o x x o o o x o o  o  o 
12 x x x x o o o o o x x  o  o 
13 x x x x o o o o o o o  x  x 
14 x x x o o o o o o o o  x  o 
15 x x x x o o o x o o o  x  x 
16 x o o o o o x x o o o  x  o 
17 x o o o o o x o o o o  x  x 
18 x o o o o x o o o o x  x  x 
19 o o o o x x o o o o x  x  o 
20 o x x x o o o o o o o  x  o 
21 x x x x x o x o o o o  x  o 
22 x x o x o x o o o o o  x  o 
23 x x o x x x o o o o o  x  x 
24 o x x x x o o o o x x  x  x 
25 x o o x x o o o o x x  x  o 
In [104]:
par(bg = 'white')
# FIXME
m1 = arima(pow, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
m1
tsdiag(m1, gof = 36)
# predict(m1, 24)
Call:
arima(x = pow, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))

Coefficients:
          ma1     sma1
      -0.4865  -0.9664
s.e.   0.0631   0.1179

sigma^2 estimated as 0.0003337:  log likelihood = 633.59,  aic = -1263.17
No description has been provided for this image
In [111]:
pow_ts = ts(pow, frequency = 12, start = c(2010, 1))
npts = 24
eotr = length(pow_ts)-npts
h = npts
freq = 12
# The order and seasonal following in combination gives a multiplicative seasonal model.
order = c(0, 1, 1)
seasonal = list(order = c(0, 1, 1), period = 12)
fixed = NULL
pow_fc_res = plot_arima_forecast_fig(
    da_ts=pow_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML',
    include.mean=T, transform.pars=NULL,
    main="Forecasts from ARIMA(0,1,1)(0,1,1)[12] for Power consumption",
    xlab="Year", ylab="Power", ylim=c(9.6, 10.1)
)
[1] 227
2027.917 ; 2032
No description has been provided for this image
In [44]:
pow_fc_tb = comb_forecast_res(pow_fc_res, pow_ts, eotr, freq)
pow_fc_tb
Forecast method: ARIMA(0,1,1)(0,1,1)[12]

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
          ma1     sma1
      -0.5018  -0.9721
s.e.   0.0658   0.1861

sigma^2 estimated as 0.0003399:  log likelihood = 569.22,  aic = -1134.43

Error measures:
                       ME       RMSE        MAE         MPE      MAPE      MASE
Training set 0.0005686413 0.01806016 0.01399088 0.005729905 0.1463576 0.4023562
                   ACF1
Training set 0.04081969

Forecasts:
         Point Forecast    Lo 80     Hi 80    Lo 95     Hi 95
Jan 2030       9.771751 9.747797  9.795706 9.735116  9.808386
Feb 2030       9.757840 9.731082  9.784598 9.716917  9.798764
Mar 2030       9.760923 9.731628  9.790219 9.716120  9.805727
Apr 2030       9.779608 9.747979  9.811238 9.731235  9.827981
May 2030       9.808923 9.775120  9.842725 9.757226  9.860619
Jun 2030       9.893825 9.857981  9.929670 9.839006  9.948645
Jul 2030       9.934974 9.897198  9.972750 9.877201  9.992748
Aug 2030       9.950545 9.910932  9.990159 9.889962 10.011129
Sep 2030       9.959941 9.918572 10.001311 9.896672 10.023211
Oct 2030       9.893905 9.850851  9.936959 9.828060  9.959750
Nov 2030       9.839844 9.795169  9.884519 9.771520  9.908168
Dec 2030       9.816089 9.769850  9.862328 9.745373  9.886806
Jan 2031       9.797054 9.748953  9.845154 9.723490  9.870617
Feb 2031       9.783142 9.733433  9.832851 9.707119  9.859165
Mar 2031       9.786225 9.734959  9.837492 9.707820  9.864631
Apr 2031       9.804911 9.752132  9.857689 9.724193  9.885629
May 2031       9.834225 9.779977  9.888473 9.751259  9.917191
Jun 2031       9.919128 9.863448  9.974807 9.833973 10.004282
Jul 2031       9.960276 9.903202 10.017351 9.872989 10.047564
Aug 2031       9.975848 9.917411 10.034284 9.886477 10.065218
Sep 2031       9.985244 9.925477 10.045011 9.893838 10.076649
Oct 2031       9.919207 9.858138  9.980276 9.825810 10.012604
Nov 2031       9.865146 9.802803  9.927490 9.769800  9.960492
Dec 2031       9.841391 9.777799  9.904984 9.744135  9.938648
A Time Series: 2 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20309.7717519.7578409.7609239.7796089.8089239.8938259.9349749.9505459.9599419.8939059.8398449.816089
20319.7970549.7831429.7862259.8049119.8342259.9191289.9602769.9758489.9852449.9192079.8651469.841391
A Time Series: 2 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
20300.018691680.020879720.022859280.024680570.026376390.027969590.029476800.030910610.032280790.033595140.034859960.03608048
20310.037533130.038788050.040003630.041183350.042330200.043446780.044535380.045598000.046636420.047652210.048646790.04962145
A Time Series: 2 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
2030 9.7696 9.7559 9.7804 9.8166 9.8245 9.8991 9.9229 9.9577 9.9685 9.9279 9.8715 9.8323
2031 9.8411 9.8128 9.8544 9.8625 9.9017 9.985710.028410.049910.0598 9.9722 9.9080 9.9060
A Time Series: 24 × 3
ForecastStd. ErrorActual
Jan 20309.7717510.01869168 9.7696
Feb 20309.7578400.02087972 9.7559
Mar 20309.7609230.02285928 9.7804
Apr 20309.7796080.02468057 9.8166
May 20309.8089230.02637639 9.8245
Jun 20309.8938250.02796959 9.8991
Jul 20309.9349740.02947680 9.9229
Aug 20309.9505450.03091061 9.9577
Sep 20309.9599410.03228079 9.9685
Oct 20309.8939050.03359514 9.9279
Nov 20309.8398440.03485996 9.8715
Dec 20309.8160890.03608048 9.8323
Jan 20319.7970540.03753313 9.8411
Feb 20319.7831420.03878805 9.8128
Mar 20319.7862250.04000363 9.8544
Apr 20319.8049110.04118335 9.8625
May 20319.8342250.04233020 9.9017
Jun 20319.9191280.04344678 9.9857
Jul 20319.9602760.0445353810.0284
Aug 20319.9758480.0455980010.0499
Sep 20319.9852440.0466364210.0598
Oct 20319.9192070.04765221 9.9722
Nov 20319.8651460.04864679 9.9080
Dec 20319.8413910.04962145 9.9060

2-7 (Equal-Weighted Index)

In [114]:
require(xts)
require(sandwich)
require(lmtest)
require(forecast)
In [13]:
find("read.table")
'package:utils'

Fit a time series model for CRSP equal-weighted index

  • Test the effect of trading days using regression model without time series error.
  • Add time series error model by using ARIMA model. The order of ARIMA can be seen from some acf and pacf plots; Also can look at eacf plots and tables. After getting a starting point for the order parameter of arima function, we can trim the number of coefs by reducing element numbers in order, until we have all coefs close to be statistically significant.
  • The logic of determining the order in arima model is as following.
    • We see that the m1$residuals have some autocorrelation in acf. We need to take diff.
    • After taking one diff, we see that the diff(m1$residuals) have some autocorrelation in both pacf and acf. We take seasonal diff.
    • After taking one seasonal diff, we see that pacf has some strong seasonality. And acf has some autocorrelation. Let's look at effect of both regular diff and seasonal diff.
    • After taking both regular and seasonal diff, we see that the autocorrelation in both pacf and acf getting stronger. We probably don't want to do this.
    • So, we need to apply seasonal diff and model the decaying autocorrelation in the pacf. That means modeling the AR process.
In [113]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = T)
da[1:5,]
A data.frame: 5 × 12
yearmomdayibmvwewspMTWRF
<int><int><int><dbl><dbl><dbl><dbl><int><int><int><int><int>
1198012-0.029126-0.020089-0.011686-0.02019600100
2198013 0.016000-0.006510-0.011628-0.00510600010
3198014-0.001969 0.013735 0.015809 0.01235500001
4198017-0.003945 0.004368 0.007013 0.00272210000
5198018 0.067327 0.019340 0.014152 0.02003601000
In [115]:
ew = da$ew * 100
ew_ts = ts(ew, frequency = 252, start = c(1980, 1, 2))
plot(
    xts(ew, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
    type = 'l', main = '', xlab = 'date', ylab = 'ew'
)
No description has been provided for this image
  • Trading day effects are all significant.
In [116]:
M = da$M
T = da$T
W = da$W
R = da$R
m1 = lm(ew ~ M + T + W + R)
summary(m1)
Call:
lm(formula = ew ~ M + T + W + R)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.2962  -0.3094   0.0533   0.3795  10.8319 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.22386    0.02155  10.389  < 2e-16 ***
M           -0.31734    0.03085 -10.286  < 2e-16 ***
T           -0.19778    0.03028  -6.532 6.94e-11 ***
W           -0.10185    0.03027  -3.365 0.000770 ***
R           -0.10294    0.03042  -3.384 0.000719 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8234 on 7314 degrees of freedom
Multiple R-squared:  0.01618,	Adjusted R-squared:  0.01564 
F-statistic: 30.06 on 4 and 7314 DF,  p-value: < 2.2e-16
  • Trading day effects are still significant. So HAC estimator doesn't change the conclusion of weekday effects.
In [17]:
find("NeweyWest")
'package:sandwich'
In [117]:
coeftest(m1, NeweyWest(m1, lag = 12, prewhite = F))
t test of coefficients:

             Estimate Std. Error  t value  Pr(>|t|)    
(Intercept)  0.223862   0.018984  11.7924 < 2.2e-16 ***
M           -0.317340   0.029024 -10.9338 < 2.2e-16 ***
T           -0.197783   0.026948  -7.3395 2.375e-13 ***
W           -0.101854   0.026331  -3.8683 0.0001106 ***
R           -0.102943   0.025381  -4.0560 5.044e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  • Residuals have serial correlation.
    • Box-Ljung test tells us the null-hypothesis can be rejected.
In [118]:
Box.test(m1$residuals, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  m1$residuals
X-squared = 514.96, df = 12, p-value < 2.2e-16
  • Residuals have some autocorrelations.
In [119]:
par(bg = 'white')
plot_time_fig(m1$residuals)
plot_time_fig(m1$residuals[1:200])
plot_time_fig(m1$residuals[(length(m1$residuals)-200):length(m1$residuals)])
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [122]:
par(mfrow = c(2, 1), bg = 'white')
pacf(m1$residuals)
acf(m1$residuals)
pacf(diff(m1$residuals))
acf(diff(m1$residuals))
pacf(diff(m1$residuals, 5))
acf(diff(m1$residuals, 5))
pacf(diff(diff(m1$residuals), 5))
acf(diff(diff(m1$residuals), 5))
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
  • Fit an ARIMA model with time series error.
In [59]:
ew_eacf_res <- perform_and_print_eacf(ew, ar.max = 25, ma.max = 12)
A data.frame: 26 × 13
0123456789101112
<I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>><I<chr>>
1 0.1991 0.0457 0.09133 0.054374 0.03882 5.56e-02 0.03096 0.02780 0.03645 0.027609 0.004854 0.06137 0.03047
2-0.0330-0.0827 0.05637-0.008304-0.01521 3.13e-02-0.01383-0.00505 0.01116 0.021850-0.000777 0.05525-0.01356
3-0.0796-0.3413 0.04200-0.010707-0.02214 2.35e-02-0.01466-0.01345 0.00726 0.001192-0.013286 0.04953-0.01297
4-0.2366-0.1846-0.18664-0.037721 0.00298-3.29e-03-0.01294 0.01206-0.00890-0.000875-0.015403 0.04869-0.00071
5-0.5000 0.0493-0.32938-0.066020-0.01221-2.66e-06-0.00280-0.00139-0.00883-0.001862-0.010022 0.01152 0.01163
6-0.4371-0.1815-0.36205-0.142458 0.02558-6.49e-03-0.00268 0.00162-0.00515-0.006756-0.014952 0.01954 0.00104
7-0.1738-0.3438-0.42397 0.155797-0.10697 4.26e-02 0.00701 0.00139-0.00703-0.004850-0.017060 0.01121 0.00826
8-0.3802-0.0588-0.29784-0.202874 0.04540 4.84e-02-0.00116 0.00208-0.01975 0.005105-0.008692 0.01378-0.01034
9-0.4783 0.1760-0.03593-0.153479 0.03060 4.32e-02 0.00780 0.06465-0.01597 0.002325-0.013141 0.00941-0.00838
10-0.4195 0.2650-0.03748-0.090770 0.14952-2.55e-01 0.03402 0.06760-0.04079-0.028994-0.007919 0.00119-0.00881
11 0.4903-0.0832-0.16157-0.186672 0.10171-2.80e-01 0.08032 0.04187-0.15953 0.021513-0.006925 0.00189 0.00126
12 0.1535-0.1568-0.33435-0.297189-0.16982-3.45e-01-0.16917-0.14890-0.14307-0.041255 0.003619 0.00257 0.00156
13-0.0354-0.2572 0.00564-0.450855 0.26370 1.88e-01-0.09841-0.17321-0.12558-0.027698 0.018525 0.00214-0.01098
14-0.1252-0.2557-0.00254-0.368854 0.13005 2.15e-01-0.19295-0.12754 0.02938-0.051806-0.007715-0.01742-0.18524
15-0.0882-0.2607-0.00773-0.253866 0.19894 2.00e-01-0.27774-0.15243-0.01056-0.109390 0.126849 0.03959-0.19531
16-0.0273-0.2508 0.01765-0.298093 0.23613 2.10e-01-0.21458-0.05457 0.09607-0.189322 0.168713 0.04839-0.20277
17 0.4257 0.4224 0.06380-0.108750 0.03701-1.99e-01 0.11547-0.13522-0.02931-0.211889-0.105417-0.12144-0.27727
18-0.5003 0.4423 0.13404-0.069572-0.01253-5.13e-02-0.01848-0.16813 0.01347-0.245807 0.013330 0.05827-0.27105
19 0.4965 0.4987 0.15593-0.196893-0.06260-5.48e-02 0.01833-0.01737-0.03233-0.240333 0.037863 0.01612-0.13779
20-0.4919 0.4729 0.13547-0.317312 0.29066-2.21e-02 0.02868-0.09510 0.00284-0.188563-0.078972-0.05367-0.09741
21 0.4248 0.3451 0.17017-0.017905-0.07052-1.03e-01 0.02755-0.03374-0.01763-0.194757-0.109101-0.04638 0.06374
22 0.3304 0.1048 0.16782-0.004319-0.05508-6.38e-02-0.04872-0.06977-0.11524-0.257878-0.064314-0.03350 0.01807
23 0.1085-0.0579 0.34436-0.000961-0.21443 2.26e-01-0.20265 0.21078 0.04499-0.238023 0.085021-0.00031-0.02534
24 0.0695-0.4345 0.34362-0.001231-0.17837-4.95e-02-0.11285 0.10954 0.08377-0.133928 0.098015 0.00314-0.09128
25-0.2373-0.0864-0.42230 0.083429-0.18267 6.07e-02-0.02557 0.19995 0.11316 0.046554 0.011300-0.04570-0.08013
26-0.4522 0.0342-0.34910-0.068288-0.15929-4.90e-02-0.04732 0.17984 0.02677 0.041395 0.008066-0.06806-0.08910
A matrix: 26 × 13 of type chr
0123456789101112
0xxxxxxxxxxoxx
1xxxooxoooooxo
2xxxooxoooooxo
3xxxxoooooooxo
4xxxxooooooooo
5xxxxxoooooooo
6xxxxxxooooooo
7xxxxxxooooooo
8xxxxxxoxooooo
9xxxxxxxxxxooo
10xxxxxxxxxoooo
11xxxxxxxxxxooo
12xxoxxxxxxxooo
13xxoxxxxxxxoox
14xxoxxxxxoxxxx
15xxoxxxxxxxxxx
16xxxxxxxxxxxxx
17xxxxoxoxoxoxx
18xxxxxxooxxxox
19xxxxxoxxoxxxx
20xxxoxxxxoxxxx
21xxxoxxxxxxxxo
22xxxoxxxxxxxox
23xxxoxxxxxxxox
24xxxxxxxxxxoxx
25xxxxxxxxxxoxx
0.0233778260111891
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12
0  x x x x x x x x x x o  x  x 
1  x x x o o x o o o o o  x  o 
2  x x x o o x o o o o o  x  o 
3  x x x x o o o o o o o  x  o 
4  x x x x o o o o o o o  o  o 
5  x x x x x o o o o o o  o  o 
6  x x x x x x o o o o o  o  o 
7  x x x x x x o o o o o  o  o 
8  x x x x x x o x o o o  o  o 
9  x x x x x x x x x x o  o  o 
10 x x x x x x x x x o o  o  o 
11 x x x x x x x x x x o  o  o 
12 x x o x x x x x x x o  o  o 
13 x x o x x x x x x x o  o  x 
14 x x o x x x x x o x x  x  x 
15 x x o x x x x x x x x  x  x 
16 x x x x x x x x x x x  x  x 
17 x x x x o x o x o x o  x  x 
18 x x x x x x o o x x x  o  x 
19 x x x x x o x x o x x  x  x 
20 x x x o x x x x o x x  x  x 
21 x x x o x x x x x x x  x  o 
22 x x x o x x x x x x x  o  x 
23 x x x o x x x x x x x  o  x 
24 x x x x x x x x x x o  x  x 
25 x x x x x x x x x x o  x  x 
In [56]:
par(bg = 'white')
m2 = arima(
    ew, order = c(2, 0, 2),
    seasonal = list(order = c(1, 0, 0), period = 5),
    xreg = da[, 8:11]
)
m2
tsdiag(m2, gof = 20)
Call:
arima(x = ew, order = c(2, 0, 2), seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11])

Coefficients:
         ar1     ar2      ma1      ma2    sar1  intercept        M        T
      0.5409  0.3211  -0.3431  -0.4035  -0.041     0.2242  -0.3171  -0.1982
s.e.  0.0950  0.0731   0.0928   0.0580   0.013     0.0239   0.0267   0.0280
            W        R
      -0.1058  -0.0987
s.e.   0.0282   0.0263

sigma^2 estimated as 0.6412:  log likelihood = -8759.03,  aic = 17538.07
No description has been provided for this image
In [58]:
par(bg = 'white')
m3 = arima(
    ew, order = c(2, 0, 3), 
    seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11]
)
m3
tsdiag(m3, gof = 20)
Call:
arima(x = ew, order = c(2, 0, 3), seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11])

Coefficients:
         ar1     ar2      ma1      ma2      ma3     sar1  intercept        M
      0.4950  0.3675  -0.2966  -0.4397  -0.0099  -0.0412     0.2242  -0.3169
s.e.  0.1635  0.1548   0.1639   0.1209   0.0296   0.0131     0.0240   0.0267
            T        W        R
      -0.1983  -0.1056  -0.0986
s.e.   0.0280   0.0282   0.0263

sigma^2 estimated as 0.6412:  log likelihood = -8758.99,  aic = 17539.97
No description has been provided for this image
In [24]:
par(bg = 'white')
m4 = arima(
    ew, order = c(4, 0, 4), 
    seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11]
)
m4
tsdiag(m4, gof = 20)
Call:
arima(x = ew, order = c(4, 0, 4), seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11])

Coefficients:
         ar1      ar2     ar3     ar4     ma1     ma2      ma3      ma4
      0.0926  -0.3079  0.8712  0.0879  0.1076  0.3207  -0.7429  -0.2188
s.e.  0.1167   0.0213  0.0516  0.0930  0.1155  0.0311   0.0638   0.0756
         sar1  intercept        M        T        W        R
      -0.0115     0.2241  -0.3161  -0.1979  -0.1066  -0.0977
s.e.   0.0138     0.0242   0.0268   0.0283   0.0285   0.0264

sigma^2 estimated as 0.6367:  log likelihood = -8733.3,  aic = 17496.59
No description has been provided for this image
In [25]:
par(bg = 'white')
m4 = arima(
    ew, order = c(3, 0, 3), 
    seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11]
)
m4
tsdiag(m4, gof = 20)
Warning message in arima(ew, order = c(3, 0, 3), seasonal = list(order = c(1, 0, :
“possible convergence problem: optim gave code = 1”
Call:
arima(x = ew, order = c(3, 0, 3), seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11])

Coefficients:
          ar1     ar2     ar3     ma1      ma2      ma3     sar1  intercept
      -0.5192  0.8231  0.3977  0.7205  -0.7018  -0.4755  -0.0350     0.2219
s.e.   0.1057  0.0431  0.0761  0.1032   0.0637   0.0585   0.0139     0.0239
            M        T        W        R
      -0.3145  -0.1955  -0.1019  -0.0952
s.e.   0.0267   0.0281   0.0283   0.0263

sigma^2 estimated as 0.64:  log likelihood = -8752.36,  aic = 17530.72
No description has been provided for this image
In [126]:
par(bg = 'white')
m5 = arima(
    ew, order = c(2, 0, 2), 
#     seasonal = list(order = c(0, 0, 1), period = 5), 
    xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
Call:
arima(x = ew, order = c(2, 0, 2), xreg = da[, 8:11])

Coefficients:
         ar1     ar2      ma1      ma2  intercept        M        T        W
      0.6226  0.2504  -0.4251  -0.3515     0.2244  -0.3181  -0.1988  -0.1053
s.e.  0.0948  0.0685   0.0931   0.0552     0.0243   0.0275   0.0291   0.0293
            R
      -0.0994
s.e.   0.0271

sigma^2 estimated as 0.6421:  log likelihood = -8764.1,  aic = 17546.19
No description has been provided for this image
In [125]:
par(bg = 'white')
m5 = arima(
    ew, order = c(2, 0, 2), 
    seasonal = list(order = c(0, 0, 1), period = 5), 
    xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
Call:
arima(x = ew, order = c(2, 0, 2), seasonal = list(order = c(0, 0, 1), period = 5), 
    xreg = da[, 8:11])

Coefficients:
         ar1     ar2      ma1      ma2     sma1  intercept        M        T
      0.5399  0.3229  -0.3421  -0.4050  -0.0426     0.2241  -0.3171  -0.1981
s.e.  0.0945  0.0734   0.0923   0.0584   0.0133     0.0239   0.0267   0.0279
            W        R
      -0.1054  -0.0983
s.e.   0.0281   0.0262

sigma^2 estimated as 0.6412:  log likelihood = -8758.85,  aic = 17537.69
No description has been provided for this image
In [128]:
par(bg = 'white')
m5 = arima(
    ew, order = c(2, 0, 2), 
    seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
Call:
arima(x = ew, order = c(2, 0, 2), seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11])

Coefficients:
         ar1     ar2      ma1      ma2    sar1  intercept        M        T
      0.5409  0.3211  -0.3431  -0.4035  -0.041     0.2242  -0.3171  -0.1982
s.e.  0.0950  0.0731   0.0928   0.0580   0.013     0.0239   0.0267   0.0280
            W        R
      -0.1058  -0.0987
s.e.   0.0282   0.0263

sigma^2 estimated as 0.6412:  log likelihood = -8759.03,  aic = 17538.07
No description has been provided for this image
In [127]:
par(bg = 'white')
m5 = arima(
    ew, order = c(2, 0, 2), 
    seasonal = list(order = c(1, 0, 1), period = 5), 
    xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
Call:
arima(x = ew, order = c(2, 0, 2), seasonal = list(order = c(1, 0, 1), period = 5), 
    xreg = da[, 8:11])

Coefficients:
         ar1     ar2     ma1      ma2    sar1     sma1  intercept        M
      0.5554  0.3206  -0.358  -0.4060  0.2901  -0.3341     0.2244  -0.3173
s.e.  0.0924  0.0738   0.090   0.0592  0.1727   0.1702     0.0239   0.0264
            T        W        R
      -0.1977  -0.1062  -0.0981
s.e.   0.0276   0.0277   0.0260

sigma^2 estimated as 0.641:  log likelihood = -8757.93,  aic = 17537.86
No description has been provided for this image

Forecast using the model

In [7]:
help(auto.arima)

Test auto.arima with a subset of time series data

In [17]:
da[6560:6570,]
A data.frame: 11 × 12
yearmomdayibmvwewspMTWRF
<int><int><int><dbl><dbl><dbl><dbl><int><int><int><int><int>
656020051223 0.003124 0.000981 0.003363 0.00042600001
656120051227-0.005870-0.009806-0.007569-0.00955301000
656220051228 0.000603 0.002857 0.003299 0.00129700100
656320051229-0.007707-0.002343 0.000072-0.00298100010
656420051230-0.002427-0.004415 0.000053-0.00488700001
65652006 1 3-0.001703 0.016428 0.010981 0.01643001000
65662006 1 4-0.001340 0.005531 0.007697 0.00367300100
65672006 1 5 0.006711-0.000332 0.002920 0.00001600010
65682006 1 6 0.029697 0.009813 0.009032 0.00939900001
65692006 1 9-0.014361 0.004303 0.007009 0.00365610000
65702006 110 0.004061 0.000978 0.004278-0.00035701000
In [19]:
sub_ew_ts = ts(ew[6565:length(ew)], frequency = 252, start = c(2006, 1, 3))
sub_da = da[6565:length(ew),]
length(sub_ew_ts)
755
In [23]:
sub_ts_fm <- auto.arima(
    sub_ew_ts,
    d = 0,
    D = 0,
    max.p = 2,
    max.q = 2,
    max.P = 1,
    max.Q = 0,
    max.order = 5,
#     nmodels = 10,
    seasonal = TRUE,
    method = "ML",
    allowmean = TRUE,
    xreg = as.matrix(sub_da[, 8:11]),
#     stepwise = FALSE,
#     parallel = TRUE,
#     num.cores = 6
)
sub_ts_fm
Series: sub_ew_ts 
Regression with ARIMA(2,0,2) errors 

Coefficients:
          ar1      ar2     ma1     ma2        M       T        W       R
      -0.9976  -0.6214  1.0519  0.5740  -0.2183  0.0709  -0.0023  -0.060
s.e.   0.2445   0.2018  0.2520  0.2289   0.1147  0.1108   0.1099   0.111

sigma^2 = 2.013:  log likelihood = -1331.43
AIC=2680.87   AICc=2681.11   BIC=2722.51
In [25]:
class(sub_ts_fm); attributes(sub_ts_fm)
  1. 'forecast_ARIMA'
  2. 'ARIMA'
  3. 'Arima'
$names
  1. 'coef'
  2. 'sigma2'
  3. 'var.coef'
  4. 'mask'
  5. 'loglik'
  6. 'aic'
  7. 'arma'
  8. 'residuals'
  9. 'call'
  10. 'series'
  11. 'code'
  12. 'n.cond'
  13. 'nobs'
  14. 'model'
  15. 'bic'
  16. 'aicc'
  17. 'xreg'
  18. 'x'
  19. 'fitted'
$class
  1. 'forecast_ARIMA'
  2. 'ARIMA'
  3. 'Arima'
In [29]:
h = 20
sub_ts_fc_res <- forecast(sub_ts_fm, h = h, xreg = as.matrix(sub_da[(dim(sub_da)[1]-h):(dim(sub_da)[1]), 8:11]))
summary(sub_ts_fc_res)
Forecast method: Regression with ARIMA(2,0,2) errors

Model Information:
Series: sub_ew_ts 
Regression with ARIMA(2,0,2) errors 

Coefficients:
          ar1      ar2     ma1     ma2        M       T        W       R
      -0.9976  -0.6214  1.0519  0.5740  -0.2183  0.0709  -0.0023  -0.060
s.e.   0.2445   0.2018  0.2520  0.2289   0.1147  0.1108   0.1099   0.111

sigma^2 = 2.013:  log likelihood = -1331.43
AIC=2680.87   AICc=2681.11   BIC=2722.51

Error measures:
                     ME     RMSE       MAE      MPE     MAPE      MASE
Training set 0.01357602 1.411296 0.8930034 106.3844 138.8498 0.7143814
                   ACF1
Training set 0.01716401

Forecasts:
         Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
2008.996   -0.116092511 -1.934400 1.702215 -2.896953 2.664768
2009.000   -0.149994315 -1.970978 1.670990 -2.934949 2.634960
2009.004    0.203478305 -1.626850 2.033807 -2.595768 3.002724
2009.008   -0.171110322 -2.005563 1.663343 -2.976664 2.634443
2009.012   -0.211340702 -2.045811 1.623129 -3.016920 2.594239
2009.016    0.170234654 -1.665515 2.005985 -2.637303 2.977772
2009.020   -0.105806788 -1.943019 1.731405 -2.915580 2.703966
2009.024   -0.018528507 -1.855994 1.818937 -2.828689 2.791632
2009.028    0.022914938 -1.814612 1.860442 -2.787340 2.833170
2009.032   -0.266953625 -2.104796 1.570888 -3.077690 2.543783
2009.036    0.105149309 -1.732857 1.943155 -2.705838 2.916137
2009.040   -0.006317078 -1.844326 1.831692 -2.817309 2.804675
2009.044   -0.077339654 -1.915387 1.760708 -2.888391 2.733712
2009.048    0.019764592 -1.818336 1.857865 -2.791368 2.830897
2009.052   -0.227265507 -2.065378 1.610847 -3.038416 2.583885
2009.056    0.067512392 -1.770601 1.905626 -2.743640 2.878665
2009.060    0.006568932 -1.831555 1.844693 -2.804599 2.817737
2009.063   -0.006801287 -1.844932 1.831329 -2.817979 2.804377
2009.067   -0.217068216 -2.055199 1.621063 -3.028247 2.594110
2009.071    0.073847161 -1.764285 1.911979 -2.737333 2.885027
2009.075   -0.006087489 -1.844221 1.832046 -2.817271 2.805096
In [31]:
par(bg = 'white')
plot(sub_ts_fc_res)
No description has been provided for this image
In [32]:
sub_ts_fc_res$mean
A Time Series:
7299
-0.116092511465873
7300
-0.149994315460839
7301
0.203478304950891
7302
-0.171110322048078
7303
-0.211340702172188
7304
0.170234654498533
7305
-0.105806788188471
7306
-0.0185285070907785
7307
0.022914938280641
7308
-0.266953625499702
7309
0.105149309488776
7310
-0.00631707820380557
7311
-0.0773396544195852
7312
0.0197645917869954
7313
-0.227265507085973
7314
0.0675123922873163
7315
0.00656893192937669
7316
-0.00680128731368499
7317
-0.217068215706301
7318
0.0738471605783947
7319
-0.00608748867555768
In [33]:
sub_ts_fc_res$lower; sub_ts_fc_res$upper
A Time Series: 21 × 2
80%95%
2008.996-1.934400-2.896953
2009.000-1.970978-2.934949
2009.004-1.626850-2.595768
2009.008-2.005563-2.976664
2009.012-2.045811-3.016920
2009.016-1.665515-2.637303
2009.020-1.943019-2.915580
2009.024-1.855994-2.828689
2009.028-1.814612-2.787340
2009.032-2.104796-3.077690
2009.036-1.732857-2.705838
2009.040-1.844326-2.817309
2009.044-1.915387-2.888391
2009.048-1.818336-2.791368
2009.052-2.065378-3.038416
2009.056-1.770601-2.743640
2009.060-1.831555-2.804599
2009.063-1.844932-2.817979
2009.067-2.055199-3.028247
2009.071-1.764285-2.737333
2009.075-1.844221-2.817271
A Time Series: 21 × 2
80%95%
2008.9961.7022152.664768
2009.0001.6709902.634960
2009.0042.0338073.002724
2009.0081.6633432.634443
2009.0121.6231292.594239
2009.0162.0059852.977772
2009.0201.7314052.703966
2009.0241.8189372.791632
2009.0281.8604422.833170
2009.0321.5708882.543783
2009.0361.9431552.916137
2009.0401.8316922.804675
2009.0441.7607082.733712
2009.0481.8578652.830897
2009.0521.6108472.583885
2009.0561.9056262.878665
2009.0601.8446932.817737
2009.0631.8313292.804377
2009.0671.6210632.594110
2009.0711.9119792.885027
2009.0751.8320462.805096

Plot forecast results from auto.arima model with a subset of time series

In [41]:
help(stopifnot)
In [112]:
npts = 20
eotr = length(sub_ew_ts)-npts
h = npts
freq = 252
xreg = as.matrix(sub_da[, 8:11])
sub_ew_fc_res = plot_auto_arima_forecast_fig(
    da_ts=sub_ew_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    xreg=xreg,
    main="Forecasts from ARIMA(1,1,0)\nfor CRSP Equal-Weighted Index", 
    xlab="Year", ylab="CRSP Equal-Weighted Index"# , ylim=c(5, 11)
)
[1] 734
2008.833 ; 2008.996
No description has been provided for this image
In [49]:
sub_ew_fc_tb = comb_forecast_res(sub_ew_fc_res, sub_ew_ts, eotr, freq)
sub_ew_fc_tb
Forecast method: Regression with ARIMA(1,1,0) errors

Model Information:
Series: tr_da_ts 
Regression with ARIMA(1,1,0) errors 

Coefficients:
          ar1    drift     xreg
      -0.4319  -0.0011  -0.1894
s.e.   0.0343   0.0444   0.1196

sigma^2 = 2.983:  log likelihood = -1441.2
AIC=2890.39   AICc=2890.45   BIC=2908.79

Error measures:
                        ME     RMSE      MAE      MPE     MAPE     MASE
Training set -0.0001707146 1.722408 1.026185 128.2297 323.7154 0.850973
                   ACF1
Training set -0.1714507

Forecasts:
         Point Forecast     Lo 80     Hi 80      Lo 95     Hi 95
2008.917   -1.298624463 -3.512011 0.9147619  -4.683707  2.086458
2008.921    0.849334056 -1.696299 3.3949670  -3.043876  4.742544
2008.925   -0.079878948 -3.124575 2.9648171  -4.736340  4.576582
2008.929    0.130539853 -3.260069 3.5211483  -5.054948  5.316028
2008.933    0.145715474 -3.590326 3.8817573  -5.568067  5.859498
2008.937    0.219419115 -3.819953 4.2587911  -5.958267  6.397105
2008.940    0.186057742 -4.140531 4.5126467  -6.430889  6.803004
2008.944    0.198936766 -4.394866 4.7927393  -6.826678  7.224551
2008.948    0.002474672 -4.844685 4.8496339  -7.410616  7.415565
2008.952    0.193378414 -4.894185 5.2809422  -7.587379  7.974136
2008.956    0.191186694 -5.126069 5.5084420  -7.940853  8.323227
2008.960    0.190603831 -5.346763 5.7279709  -8.278068  8.659276
2008.964    0.189326118 -5.559757 5.9384089  -8.603137  8.981789
2008.968   -0.001021811 -5.954285 5.9522417  -9.105753  9.103709
2008.972    0.187241279 -5.963433 6.3379156  -9.219403  9.593886
2008.976    0.186190032 -6.155751 6.5281314  -9.512972  9.885352
2008.980    0.185114609 -6.342493 6.7127217  -9.797999 10.168228
2008.984   -0.005320687 -6.713456 6.7028151 -10.264529 10.253888
2008.988    0.182980137 -6.700952 7.0669120 -10.345085 10.711045
2008.992    0.181912593 -6.873436 7.2372616 -10.608313 10.972138
A Time Series:
  1. -1.29862446268606
  2. 0.849334055523345
  3. -0.0798789475238652
  4. 0.13053985253886
  5. 0.145715473588767
  6. 0.21941911505262
  7. 0.186057742363957
  8. 0.198936765953585
  9. 0.00247467231699205
  10. 0.193378414374947
  11. 0.19118669434306
  12. 0.190603831226609
  13. 0.18932611758486
  14. -0.0010218112634669
  15. 0.187241279177729
  16. 0.186190032326165
  17. 0.18511460931776
  18. -0.0053206869971778
  19. 0.182980136650425
  20. 0.181912593173573
A Time Series:
  1. 1.72711452400449
  2. 1.98636795601448
  3. 2.37578891756174
  4. 2.64570582602999
  5. 2.91524892767082
  6. 3.15193872378606
  7. 3.37605527893742
  8. 3.5845631759321
  9. 3.7822584655444
  10. 3.96984714598476
  11. 4.14907632789653
  12. 4.32083048132734
  13. 4.48603312615497
  14. 4.64535623790657
  15. 4.79939669948165
  16. 4.94864311134307
  17. 5.09351889303796
  18. 5.23438612895216
  19. 5.37156054784085
  20. 5.50531804706904
A Time Series:
  1. 1.496
  2. -2.2343
  3. 2.695
  4. 3.8066
  5. -1.5286
  6. 2.2634
  7. -2.7784
  8. 1.9321
  9. -2.1021
  10. 4.7713
  11. 1.1627
  12. -1.1401
  13. 0.3764
  14. -1.2094
  15. -0.8332
  16. 0.5254
  17. 1.1629
  18. -1.6514
  19. 2.1692
  20. 3.6731
A Time Series: 20 × 3
ForecastStd. ErrorActual
2008.917-1.2986244631.727115 1.4960
2008.921 0.8493340561.986368-2.2343
2008.925-0.0798789482.375789 2.6950
2008.929 0.1305398532.645706 3.8066
2008.933 0.1457154742.915249-1.5286
2008.937 0.2194191153.151939 2.2634
2008.940 0.1860577423.376055-2.7784
2008.944 0.1989367663.584563 1.9321
2008.948 0.0024746723.782258-2.1021
2008.952 0.1933784143.969847 4.7713
2008.956 0.1911866944.149076 1.1627
2008.960 0.1906038314.320830-1.1401
2008.964 0.1893261184.486033 0.3764
2008.968-0.0010218114.645356-1.2094
2008.972 0.1872412794.799397-0.8332
2008.976 0.1861900324.948643 0.5254
2008.980 0.1851146095.093519 1.1629
2008.984-0.0053206875.234386-1.6514
2008.988 0.1829801375.371561 2.1692
2008.992 0.1819125935.505318 3.6731
In [62]:
class(sub_ew_fc_res); attributes(sub_ew_fc_res);
class(sub_ew_fc_res$model); attributes(sub_ew_fc_res$model);
'forecast'
$names
  1. 'method'
  2. 'model'
  3. 'level'
  4. 'mean'
  5. 'lower'
  6. 'upper'
  7. 'x'
  8. 'series'
  9. 'fitted'
  10. 'residuals'
$class
'forecast'
  1. 'forecast_ARIMA'
  2. 'ARIMA'
  3. 'Arima'
$names
  1. 'coef'
  2. 'sigma2'
  3. 'var.coef'
  4. 'mask'
  5. 'loglik'
  6. 'aic'
  7. 'arma'
  8. 'residuals'
  9. 'call'
  10. 'series'
  11. 'code'
  12. 'n.cond'
  13. 'nobs'
  14. 'model'
  15. 'xreg'
  16. 'bic'
  17. 'aicc'
  18. 'x'
  19. 'fitted'
$class
  1. 'forecast_ARIMA'
  2. 'ARIMA'
  3. 'Arima'
In [59]:
sub_ew_fc_res$model
Series: tr_da_ts 
Regression with ARIMA(1,1,0) errors 

Coefficients:
          ar1    drift     xreg
      -0.4319  -0.0011  -0.1894
s.e.   0.0343   0.0444   0.1196

sigma^2 = 2.983:  log likelihood = -1441.2
AIC=2890.39   AICc=2890.45   BIC=2908.79

All data takes too long time to run auto.arima

Dummy example for auto.arima

In [50]:
# Assuming 'forecast' package is installed and loaded
library(forecast)

# Simulate some time series data
set.seed(123)
ts_data <- ts(rnorm(100), frequency = 12, start = c(2000, 1))

# Simulate exogenous regressor data for model fitting
xreg_train <- matrix(rnorm(200), ncol = 2) # Two exogenous regressors

# Fit an ARIMA model with exogenous regressors
# fit <- auto.arima(ts_data, xreg = xreg_train)
fit <- auto.arima(ts_data, xreg = xreg_train)

# Simulate future values of the exogenous regressors for forecasting
xreg_future <- matrix(rnorm(24), ncol = 2) # Two exogenous regressors for 12 future periods

# Forecast using the model and future values of exogenous regressors
forecast_result <- forecast(fit, xreg = xreg_future)

# Plot the forecast
par(bg = 'white')
plot(forecast_result)
No description has been provided for this image
In [57]:
npts = 24
eotr = length(ts_data)-npts
h = npts
freq = 12
xreg = xreg_train
sim_data_fc_res = plot_auto_arima_forecast_fig(
    da_ts=ts_data, eotr=eotr, h=h, npts=npts, frequency=freq, 
    xreg=xreg,
    main="Forecasts from ARIMA(0,0,0) for Simulated Data", 
    xlab="Year", ylab="Simulated Data"# , ylim=c(5, 11)
)
[1] 76
2004.25 ; 2008.333
No description has been provided for this image

2-8 (S&P)

In [141]:
require(xts)
In [145]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = TRUE)
dim(da); da[1:5,]
  1. 7319
  2. 12
A data.frame: 5 × 12
yearmomdayibmvwewspMTWRF
<int><int><int><dbl><dbl><dbl><dbl><int><int><int><int><int>
1198012-0.029126-0.020089-0.011686-0.02019600100
2198013 0.016000-0.006510-0.011628-0.00510600010
3198014-0.001969 0.013735 0.015809 0.01235500001
4198017-0.003945 0.004368 0.007013 0.00272210000
5198018 0.067327 0.019340 0.014152 0.02003601000
In [155]:
plot_pacf_acf(sp, lag.max = 20, main = "SP")
plot_pacf_acf(diff(sp), lag.max = 20, main = "diff(SP)")
plot_pacf_acf(diff(sp, 5), lag.max = 20, main = "diff(SP, 5)")
plot_pacf_acf(diff(diff(sp), 5), lag.max = 20, main = "diff(diff(SP), 5)")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

(a)

  • No Friday effect.
In [159]:
sp = da$sp * 100
sp_ts = ts(sp, frequency = 252, start = c(1980, 1, 2))
plot(
    xts(sp, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
    type = 'l', main = '', xlab = 'date', ylab = 'sp'
)
M = da$M
T = da$T
W = da$W
R = da$R
F = da$F
No description has been provided for this image
In [162]:
sp_m1 = lm(sp ~ M + T + W + R + F + 0)
summary(sp_m1);
Call:
lm(formula = sp ~ M + T + W + R + F + 0)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.4627  -0.5214   0.0142   0.5379  11.5842 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)  
M -0.004214   0.029944  -0.141   0.8881  
T  0.074036   0.028855   2.566   0.0103 *
W  0.063973   0.028836   2.219   0.0265 *
R  0.006432   0.029128   0.221   0.8252  
F  0.032896   0.029228   1.125   0.2604  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.117 on 7314 degrees of freedom
Multiple R-squared:  0.001753,	Adjusted R-squared:  0.00107 
F-statistic: 2.568 on 5 and 7314 DF,  p-value: 0.02501
In [163]:
Box.test(sp_m1$residuals, lag = 12, type = 'Ljung');
	Box-Ljung test

data:  sp_m1$residuals
X-squared = 67.025, df = 12, p-value = 1.149e-09
In [220]:
sp_m1 = lm(sp ~ 1 + T + W)
summary(sp_m1);
Call:
lm(formula = sp ~ 1 + T + W)

Residuals:
     Min       1Q   Median       3Q      Max 
-20.4788  -0.5215   0.0150   0.5353  11.5681 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.01195    0.01699   0.703   0.4819  
T            0.06209    0.03348   1.854   0.0637 .
W            0.05203    0.03347   1.555   0.1201  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.117 on 7316 degrees of freedom
Multiple R-squared:  0.0006394,	Adjusted R-squared:  0.0003662 
F-statistic: 2.341 on 2 and 7316 DF,  p-value: 0.09635

(b)

  • There are some serial correlation in the residuals.
In [164]:
plot_pacf_acf(sp_m1$residuals, lag.max = 20, main = 'obs')
plot_pacf_acf(diff(sp_m1$residuals), lag.max = 20, main = 'diff(obs)')
plot_pacf_acf(diff(sp_m1$residuals, 5), lag.max = 20, main = "diff(obs, 5)")
plot_pacf_acf(diff(diff(sp_m1$residuals), 5), lag.max = 20, main = "diff(diff(obs), 5)")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [169]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(0, 0, 1),
#     seasonal = list(order = c(1, 0, 0), period = 5),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Call:
arima(x = sp, order = c(0, 0, 1), xreg = da[, 8:11])

Coefficients:
          ma1  intercept        M       T       W        R
      -0.0191     0.0329  -0.0371  0.0411  0.0311  -0.0265
s.e.   0.0125     0.0292   0.0422  0.0411  0.0411   0.0416

sigma^2 estimated as 1.246:  log likelihood = -11190.07,  aic = 22392.13
No description has been provided for this image
In [170]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(1, 0, 0),
#     seasonal = list(order = c(1, 0, 0), period = 5),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Call:
arima(x = sp, order = c(1, 0, 0), xreg = da[, 8:11])

Coefficients:
          ar1  intercept        M       T       W        R
      -0.0167     0.0329  -0.0371  0.0411  0.0312  -0.0269
s.e.   0.0117     0.0292   0.0422  0.0411  0.0410   0.0416

sigma^2 estimated as 1.246:  log likelihood = -11190.21,  aic = 22392.42
No description has been provided for this image
In [171]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(2, 0, 0),
#     seasonal = list(order = c(1, 0, 0), period = 5),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Call:
arima(x = sp, order = c(2, 0, 0), xreg = da[, 8:11])

Coefficients:
          ar1      ar2  intercept        M       T       W        R
      -0.0178  -0.0627     0.0320  -0.0363  0.0426  0.0305  -0.0238
s.e.   0.0117   0.0117     0.0292   0.0421  0.0422  0.0423   0.0415

sigma^2 estimated as 1.241:  log likelihood = -11175.83,  aic = 22365.66
No description has been provided for this image
In [172]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(2, 0, 0),
    seasonal = list(order = c(1, 0, 0), period = 5),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Call:
arima(x = sp, order = c(2, 0, 0), seasonal = list(order = c(1, 0, 0), period = 5), 
    xreg = da[, 8:11])

Coefficients:
          ar1      ar2     sar1  intercept        M       T       W        R
      -0.0182  -0.0630  -0.0120     0.0320  -0.0363  0.0428  0.0303  -0.0234
s.e.   0.0117   0.0117   0.0117     0.0289   0.0417  0.0418  0.0419   0.0411

sigma^2 estimated as 1.241:  log likelihood = -11175.31,  aic = 22366.61
No description has been provided for this image
In [174]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(3, 0, 1),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Call:
arima(x = sp, order = c(3, 0, 1), xreg = da[, 8:11])

Coefficients:
          ar1      ar2      ar3     ma1  intercept        M       T       W
      -0.4878  -0.0708  -0.0167  0.4708     0.0321  -0.0363  0.0426  0.0302
s.e.   0.2230   0.0136   0.0196  0.2225     0.0292   0.0422  0.0420  0.0421
            R
      -0.0238
s.e.   0.0416

sigma^2 estimated as 1.241:  log likelihood = -11175.06,  aic = 22368.12
No description has been provided for this image
In [195]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(4, 0, 0),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Call:
arima(x = sp, order = c(4, 0, 0), xreg = da[, 8:11])

Coefficients:
          ar1      ar2     ar3      ar4  intercept        M       T       W
      -0.0176  -0.0651  0.0016  -0.0399     0.0322  -0.0367  0.0426  0.0300
s.e.   0.0117   0.0117  0.0117   0.0117     0.0291   0.0425  0.0419  0.0421
            R
      -0.0243
s.e.   0.0419

sigma^2 estimated as 1.239:  log likelihood = -11169.97,  aic = 22357.95
No description has been provided for this image
In [196]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(4, 0, 0),
    fixed = c(NA, NA, 0, NA, 0, 0, 0, 0, 0),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = sp, order = c(4, 0, 0), xreg = da[, 8:11], fixed = c(NA, NA, 0, NA, 
    0, 0, 0, 0, 0))

Coefficients:
          ar1      ar2  ar3      ar4  intercept  M  T  W  R
      -0.0166  -0.0643    0  -0.0388          0  0  0  0  0
s.e.   0.0117   0.0117    0   0.0117          0  0  0  0  0

sigma^2 estimated as 1.242:  log likelihood = -11177.18,  aic = 22360.37
No description has been provided for this image
In [198]:
par(bg = 'white')
sp_m2 = arima(
    sp, order = c(5, 0, 0),
    fixed = c(NA, NA, 0, NA, 0, 0, 0, 0, 0, 0),
    xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = sp, order = c(5, 0, 0), xreg = da[, 8:11], fixed = c(NA, NA, 0, NA, 
    0, 0, 0, 0, 0, 0))

Coefficients:
          ar1      ar2  ar3      ar4  ar5  intercept  M  T  W  R
      -0.0166  -0.0643    0  -0.0388    0          0  0  0  0  0
s.e.   0.0117   0.0117    0   0.0117    0          0  0  0  0  0

sigma^2 estimated as 1.242:  log likelihood = -11177.18,  aic = 22360.37
No description has been provided for this image
In [ ]:

In [200]:
par(bg = 'white')
sp_m3 = arima(
    sp, order = c(2, 0, 0),
    xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
Call:
arima(x = sp, order = c(2, 0, 0), xreg = da[, 8:11])

Coefficients:
          ar1      ar2  intercept        M       T       W        R
      -0.0178  -0.0627     0.0320  -0.0363  0.0426  0.0305  -0.0238
s.e.   0.0117   0.0117     0.0292   0.0421  0.0422  0.0423   0.0415

sigma^2 estimated as 1.241:  log likelihood = -11175.83,  aic = 22365.66
No description has been provided for this image
In [201]:
par(bg = 'white')
sp_m3 = arima(
    sp, order = c(2, 0, 2),
    xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
Call:
arima(x = sp, order = c(2, 0, 2), xreg = da[, 8:11])

Coefficients:
         ar1     ar2      ma1      ma2  intercept        M       T       W
      0.0737  0.3057  -0.0899  -0.3717     0.0322  -0.0361  0.0426  0.0304
s.e.  0.1258  0.1043   0.1225   0.1014     0.0290   0.0423  0.0423  0.0424
            R
      -0.0242
s.e.   0.0417

sigma^2 estimated as 1.24:  log likelihood = -11171.52,  aic = 22361.04
No description has been provided for this image
In [203]:
par(bg = 'white')
sp_m3 = arima(
    sp, order = c(2, 0, 2),
    fixed = c(0, NA, 0, NA, 0, 0, 0, 0, 0),
    xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = sp, order = c(2, 0, 2), xreg = da[, 8:11], fixed = c(0, NA, 0, NA, 
    0, 0, 0, 0, 0))

Coefficients:
      ar1     ar2  ma1      ma2  intercept  M  T  W  R
        0  0.2655    0  -0.3324          0  0  0  0  0
s.e.    0  0.1054    0   0.1028          0  0  0  0  0

sigma^2 estimated as 1.243:  log likelihood = -11180.13,  aic = 22364.27
No description has been provided for this image
In [ ]:

In [ ]:

In [213]:
sp_m4 = arima(sp, order = c(2, 0, 2), xreg = da[, 8:12], include.mean = FALSE)
sp_m4;
par(bg = 'white')
tsdiag(sp_m4, gof = 20)
Call:
arima(x = sp, order = c(2, 0, 2), xreg = da[, 8:12], include.mean = FALSE)

Coefficients:
         ar1     ar2      ma1      ma2        M       T       W       R       F
      0.0737  0.3057  -0.0899  -0.3717  -0.0039  0.0748  0.0626  0.0080  0.0322
s.e.  0.1258  0.1043   0.1225   0.1014   0.0297  0.0286  0.0286  0.0289  0.0290

sigma^2 estimated as 1.24:  log likelihood = -11171.52,  aic = 22361.04
No description has been provided for this image
  • Adding T and W, b/c they show significance in linear model.
In [215]:
sp_m4 = arima(sp, order = c(2, 0, 2), xreg = da[, c(9, 10, 12)], include.mean = FALSE,
    fixed = c(0, NA, 0, NA, NA, NA, NA))
sp_m4;
par(bg = 'white')
tsdiag(sp_m4, gof = 20)
Box.test(sp_m4$residuals, lag = 10, type = 'Ljung')
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = sp, order = c(2, 0, 2), xreg = da[, c(9, 10, 12)], include.mean = FALSE, 
    fixed = c(0, NA, 0, NA, NA, NA, NA))

Coefficients:
      ar1     ar2  ma1      ma2       T       W       F
        0  0.2736    0  -0.3413  0.0751  0.0619  0.0321
s.e.    0  0.1041    0   0.1014  0.0287  0.0286  0.0291

sigma^2 estimated as 1.24:  log likelihood = -11173.17,  aic = 22356.34
	Box-Ljung test

data:  sp_m4$residuals
X-squared = 16.368, df = 10, p-value = 0.08957
No description has been provided for this image
In [ ]:

2-9 (IBM)

In [217]:
require(xts)
In [219]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = TRUE)
da[1:5,]
A data.frame: 5 × 12
yearmomdayibmvwewspMTWRF
<int><int><int><dbl><dbl><dbl><dbl><int><int><int><int><int>
1198012-0.029126-0.020089-0.011686-0.02019600100
2198013 0.016000-0.006510-0.011628-0.00510600010
3198014-0.001969 0.013735 0.015809 0.01235500001
4198017-0.003945 0.004368 0.007013 0.00272210000
5198018 0.067327 0.019340 0.014152 0.02003601000

(a)

  • No Friday effect.
In [221]:
ibm = da$ibm * 100
plot(xts(ibm, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
    type = 'l', main = '', xlab = 'date', ylab = 'ibm')
M = da$M
T = da$T
W = da$W
R = da$R
F = da$F
No description has been provided for this image
In [223]:
ibm_m1 = lm(ibm ~ M + T + W + R + F + 0)
summary(ibm_m1)
Call:
lm(formula = ibm ~ M + T + W + R + F + 0)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.1629  -0.9290  -0.0036   0.8840  13.1619 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
M  0.199941   0.047850   4.178 2.97e-05 ***
T  0.143942   0.046109   3.122   0.0018 ** 
W -0.036134   0.046079  -0.784   0.4330    
R  0.001708   0.046546   0.037   0.9707    
F -0.059021   0.046706  -1.264   0.2064    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.785 on 7314 degrees of freedom
Multiple R-squared:  0.004006,	Adjusted R-squared:  0.003325 
F-statistic: 5.884 on 5 and 7314 DF,  p-value: 1.966e-05
In [225]:
Box.test(ibm_m1$residuals, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  ibm_m1$residuals
X-squared = 16.845, df = 12, p-value = 0.1555

(b)

  • There is no much serial correlation in the residuals.
In [231]:
lag.max = 20
par(mfrow = c(2, 1), bg = 'white')
pacf(ibm_m1$residuals, lag.max =lag.max)
acf(ibm_m1$residuals, lag.max =lag.max)
pacf(diff(ibm_m1$residuals), lag.max =lag.max)
acf(diff(ibm_m1$residuals), lag.max =lag.max)
pacf(diff(ibm_m1$residuals, 5), lag.max =lag.max)
acf(diff(ibm_m1$residuals, 5), lag.max =lag.max)
pacf(diff(diff(ibm_m1$residuals), 5), lag.max =lag.max)
acf(diff(diff(ibm_m1$residuals), 5), lag.max =lag.max)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

(c)

In [228]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, 8:12], include.mean = FALSE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
Call:
arima(x = ibm, order = c(0, 0, 1), xreg = da[, 8:12], include.mean = FALSE)

Coefficients:
          ma1       M       T        W       R        F
      -0.0321  0.2000  0.1439  -0.0361  0.0017  -0.0590
s.e.   0.0118  0.0478  0.0461   0.0461  0.0465   0.0467

sigma^2 estimated as 3.179:  log likelihood = -14618.23,  aic = 29248.46
No description has been provided for this image
In [229]:
Box.test(ibm_m2$residuals, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  ibm_m2$residuals
X-squared = 9.4174, df = 12, p-value = 0.6669
In [234]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = TRUE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
Call:
arima(x = ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = TRUE)

Coefficients:
          ma1  intercept      M       T
      -0.0323    -0.0311  0.231  0.1750
s.e.   0.0118     0.0262  0.055  0.0535

sigma^2 estimated as 3.18:  log likelihood = -14618.63,  aic = 29245.27
No description has been provided for this image
In [233]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = FALSE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
Call:
arima(x = ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = FALSE)

Coefficients:
          ma1       M       T
      -0.0322  0.1999  0.1439
s.e.   0.0118  0.0478  0.0461

sigma^2 estimated as 3.18:  log likelihood = -14619.33,  aic = 29244.65
No description has been provided for this image

2-10

  • pp10: Test statistics for skewness and kurtosis.
  • They are all significant.
In [237]:
require(fBasics)
Loading required package: fBasics


Attaching package: ‘fBasics’


The following objects are masked from ‘package:TSA’:

    kurtosis, skewness


In [235]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da2 = read.table("../AFTS_sol/data/w-Baa.txt", header = FALSE)
da1[1:5,]; da2[1:5,];
A data.frame: 5 × 4
V1V2V3V4
<int><int><int><dbl>
119621 54.43
219621124.42
319621194.42
419621264.41
519622 24.42
A data.frame: 5 × 4
V1V2V3V4
<int><int><int><dbl>
119621 55.11
219621125.09
319621195.08
419621265.08
519622 25.07
In [238]:
Aaa = da1[, 4]
Baa = da2[, 4]
apply(cbind(Aaa, Baa), 2, basicStats)
ts = skewness(Aaa) / sqrt(6 / length(Aaa))
cat("ts =", ts, "\n")
ps = (1 - pnorm(ts)) * 2
cat("ps =", ps, "\n")
tk = kurtosis(Aaa) / sqrt(24 / length(Aaa))
cat("tk =", tk, "\n")
pk = (1 - pnorm(tk)) * 2
cat("pk =", pk, "\n")
ts = skewness(Baa) / sqrt(6 / length(Baa))
cat("ts =", ts, "\n")
ps = (1 - pnorm(ts)) * 2
cat("ps =", ps, "\n")
tk = kurtosis(Baa) / sqrt(24 / length(Baa))
cat("tk =", tk, "\n")
pk = (1 - pnorm(tk)) * 2
cat("pk =", pk, "\n")
$Aaa
A data.frame: 16 × 1
X..newX..i
<dbl>
nobs 2467.000000
NAs 0.000000
Minimum 4.190000
Maximum 15.850000
1. Quartile 5.985000
3. Quartile 8.930000
Mean 7.830109
Median 7.540000
Sum19316.880000
SE Mean 0.048697
LCL Mean 7.734618
UCL Mean 7.925601
Variance 5.850323
Stdev 2.418744
Skewness 0.857092
Kurtosis 0.578605
$Baa
A data.frame: 16 × 1
X..newX..i
<dbl>
nobs 2467.000000
NAs 0.000000
Minimum 4.780000
Maximum 17.290000
1. Quartile 6.990000
3. Quartile 10.200000
Mean 8.847122
Median 8.350000
Sum21825.850000
SE Mean 0.054704
LCL Mean 8.739852
UCL Mean 8.954392
Variance 7.382486
Stdev 2.717073
Skewness 0.929779
Kurtosis 0.760896
ts = 17.37946 
ps = 0 
tk = 5.866262 
pk = 4.457306e-09 
ts = 18.85335 
ps = 0 
tk = 7.714438 
pk = 1.221245e-14 

2-11

In [239]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da1[1:5,];
A data.frame: 5 × 4
V1V2V3V4
<int><int><int><dbl>
119621 54.43
219621124.42
319621194.42
419621264.41
519622 24.42
In [252]:
aaa = da1[,4]
aaa_ts = ts(aaa, frequency = 52, start = c(1962, 1))
lg_aaa = log(aaa)
lg_aaa_ts = ts(lg_aaa, frequency = 52, start = c(1962, 1))

Fitting solution

In [253]:
sub_sz = 100
plot_time_fig(lg_aaa_ts, main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
plot_time_fig(lg_aaa_ts[1:sub_sz], main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
plot_time_fig(lg_aaa_ts[(length(lg_aaa_ts)-sub_sz):length(lg_aaa_ts)], main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [243]:
sub_sz = 100
plot_time_fig(aaa_ts, main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
plot_time_fig(aaa_ts[1:sub_sz], main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
plot_time_fig(aaa_ts[(length(aaa_ts)-sub_sz):length(aaa_ts)], main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [250]:
lag.max = 52
par(mfrow = c(2, 1), bg = 'white')
pacf(aaa, lag.max = lag.max)
acf(aaa, lag.max = lag.max)
pacf(diff(aaa), lag.max = lag.max)
acf(diff(aaa), lag.max = lag.max)
pacf(diff(aaa, 4), lag.max = lag.max)
acf(diff(aaa, 4), lag.max = lag.max)
pacf(diff(aaa, 13), lag.max = lag.max)
acf(diff(aaa, 13), lag.max = lag.max)
pacf(diff(aaa, 26), lag.max = lag.max)
acf(diff(aaa, 26), lag.max = lag.max)
pacf(diff(diff(aaa), 26), lag.max = lag.max)
acf(diff(diff(aaa), 26), lag.max = lag.max)
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [ ]:

In [256]:
par(bg = 'white')
aaa_m1 = arima(
    aaa, order = c(0, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m1
tsdiag(aaa_m1, gof = 52)
Call:
arima(x = aaa, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26))

Coefficients:
         ma1     sma1
      0.3674  -0.9807
s.e.  0.0182   0.0111

sigma^2 estimated as 0.00823:  log likelihood = 2351.57,  aic = -4699.13
No description has been provided for this image
In [259]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(1, 1, 0), 
    seasonal = list(order = c(1, 1, 0), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Call:
arima(x = aaa, order = c(1, 1, 0), seasonal = list(order = c(1, 1, 0), period = 26))

Coefficients:
         ar1     sar1
      0.3826  -0.5206
s.e.  0.0187   0.0174

sigma^2 estimated as 0.01159:  log likelihood = 1971.92,  aic = -3939.84
No description has been provided for this image
In [258]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(1, 1, 1), 
    seasonal = list(order = c(1, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Call:
arima(x = aaa, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 26))

Coefficients:
         ar1     ma1     sar1     sma1
      0.2695  0.1155  -0.0004  -0.9803
s.e.  0.0660  0.0701   0.0216   0.0118

sigma^2 estimated as 0.008178:  log likelihood = 2359.52,  aic = -4711.04
No description has been provided for this image
In [260]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(1, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Call:
arima(x = aaa, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26))

Coefficients:
         ar1     ma1     sma1
      0.2694  0.1156  -0.9804
s.e.  0.0660  0.0701   0.0114

sigma^2 estimated as 0.008178:  log likelihood = 2359.52,  aic = -4713.04
No description has been provided for this image
In [261]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(2, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Call:
arima(x = aaa, order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26))

Coefficients:
          ar1     ar2     ma1     sma1
      -0.2492  0.1869  0.6398  -0.9798
s.e.   0.1162  0.0514  0.1127   0.0112

sigma^2 estimated as 0.008156:  log likelihood = 2363.16,  aic = -4718.31
No description has been provided for this image
In [262]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(3, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Call:
arima(x = aaa, order = c(3, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26))

Coefficients:
         ar1      ar2     ar3      ma1     sma1
      0.5016  -0.1024  0.0872  -0.1211  -0.9819
s.e.  0.2284   0.0884  0.0205   0.2292   0.0117

sigma^2 estimated as 0.008117:  log likelihood = 2367.75,  aic = -4725.5
No description has been provided for this image
In [263]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(4, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Warning message in sqrt(diag(x$var.coef)):
“NaNs produced”
Call:
arima(x = aaa, order = c(4, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26))

Coefficients:
         ar1     ar2     ar3     ar4     ma1     sma1
      0.1604  0.0277  0.0713  0.0197  0.2204  -0.9820
s.e.     NaN     NaN  0.0038     NaN     NaN   0.0118

sigma^2 estimated as 0.008117:  log likelihood = 2367.61,  aic = -4723.22
No description has been provided for this image
In [265]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(9, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Call:
arima(x = aaa, order = c(9, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26))

Coefficients:
         ar1      ar2     ar3      ar4     ar5      ar6      ar7     ar8
      0.6287  -0.1528  0.0978  -0.0187  0.0431  -0.0462  -0.0413  0.0404
s.e.  0.3546   0.1361  0.0312   0.0392  0.0242   0.0288   0.0282  0.0282
          ar9      ma1     sma1
      -0.0611  -0.2502  -0.9809
s.e.   0.0205   0.3556   0.0114

sigma^2 estimated as 0.00804:  log likelihood = 2379.85,  aic = -4737.71
No description has been provided for this image
In [266]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(9, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26),
    fixed = c(NA, 0, NA, 0, NA, NA, NA, NA, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = aaa, order = c(9, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26), 
    fixed = c(NA, 0, NA, 0, NA, NA, NA, NA, NA, 0, NA))

Coefficients:
         ar1  ar2    ar3  ar4     ar5      ar6      ar7     ar8      ar9  ma1
      0.3604    0  0.066    0  0.0413  -0.0353  -0.0539  0.0296  -0.0570    0
s.e.  0.0189    0  0.019    0  0.0203   0.0216   0.0216  0.0216   0.0202    0
         sma1
      -0.9805
s.e.   0.0115

sigma^2 estimated as 0.008068:  log likelihood = 2375.89,  aic = -4735.79
No description has been provided for this image
In [267]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(9, 1, 1), 
    seasonal = list(order = c(0, 1, 1), period = 26),
    fixed = c(NA, 0, NA, 0, NA, 0, NA, 0, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = aaa, order = c(9, 1, 1), seasonal = list(order = c(0, 1, 1), period = 26), 
    fixed = c(NA, 0, NA, 0, NA, 0, NA, 0, NA, 0, NA))

Coefficients:
         ar1  ar2     ar3  ar4     ar5  ar6      ar7  ar8      ar9  ma1
      0.3577    0  0.0656    0  0.0313    0  -0.0559    0  -0.0490    0
s.e.  0.0189    0  0.0190    0  0.0190    0   0.0190    0   0.0189    0
         sma1
      -0.9807
s.e.   0.0115

sigma^2 estimated as 0.008084:  log likelihood = 2373.45,  aic = -4734.91
No description has been provided for this image
In [273]:
par(bg = 'white')
aaa_m2 = arima(aaa, order = c(9, 1, 1))
aaa_m2
tsdiag(aaa_m2, gof = 52)
Call:
arima(x = aaa, order = c(9, 1, 1))

Coefficients:
         ar1      ar2     ar3      ar4     ar5      ar6      ar7     ar8
      0.6627  -0.1655  0.1000  -0.0194  0.0384  -0.0426  -0.0427  0.0415
s.e.  0.3576   0.1364  0.0314   0.0396  0.0244   0.0284   0.0278  0.0286
          ar9      ma1
      -0.0601  -0.2866
s.e.   0.0206   0.3586

sigma^2 estimated as 0.007928:  log likelihood = 2465.16,  aic = -4910.32
No description has been provided for this image
In [276]:
par(bg = 'white')
aaa_m2 = arima(
    aaa, order = c(9, 1, 0),
    fixed = c(NA, 0, NA, 0, 0, 0, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Call:
arima(x = aaa, order = c(9, 1, 0), fixed = c(NA, 0, NA, 0, 0, 0, NA, 0, NA))

Coefficients:
         ar1  ar2     ar3  ar4  ar5  ar6      ar7  ar8      ar9
      0.3580    0  0.0700    0    0    0  -0.0532    0  -0.0466
s.e.  0.0187    0  0.0188    0    0    0   0.0188    0   0.0187

sigma^2 estimated as 0.007976:  log likelihood = 2457.81,  aic = -4907.63
No description has been provided for this image
In [ ]:

Alternative solution

In [268]:
plot(
    xts(aaa, order.by = as.Date(paste(da1[, 1], da1[, 2], da1[, 3], sep = '-'))),
    type = 'l', main = '', xlab = 'date', ylab = 'Aaa'
)
No description has been provided for this image
In [270]:
cat("order =", ar(diff(aaa), method = 'mle')$order, "\n")
order = 9 
In [271]:
aaa_al_m1 = arima(aaa, order = c(9, 1, 0))
aaa_al_m1;
par(bg = 'white')
tsdiag(aaa_al_m1, gof = 52)
Call:
arima(x = aaa, order = c(9, 1, 0))

Coefficients:
         ar1      ar2     ar3     ar4     ar5      ar6      ar7     ar8
      0.3772  -0.0579  0.0843  0.0054  0.0392  -0.0314  -0.0531  0.0274
s.e.  0.0201   0.0215  0.0215  0.0216  0.0215   0.0215   0.0215  0.0215
          ar9
      -0.0592
s.e.   0.0201

sigma^2 estimated as 0.00793:  log likelihood = 2464.86,  aic = -4911.72
No description has been provided for this image
In [291]:
aaa_al_m1 = arima(
    aaa, order = c(9, 1, 0),
    seasonal = list(order = c(1, 0, 0), period = 13)
)
aaa_al_m1;
par(bg = 'white')
tsdiag(aaa_al_m1, gof = 52)
Call:
arima(x = aaa, order = c(9, 1, 0), seasonal = list(order = c(1, 0, 0), period = 13))

Coefficients:
         ar1      ar2     ar3     ar4     ar5      ar6      ar7     ar8
      0.3771  -0.0580  0.0844  0.0053  0.0393  -0.0316  -0.0529  0.0275
s.e.  0.0201   0.0215  0.0215  0.0216  0.0215   0.0216   0.0215  0.0215
          ar9     sar1
      -0.0589  -0.0024
s.e.   0.0202   0.0205

sigma^2 estimated as 0.00793:  log likelihood = 2464.87,  aic = -4909.73
No description has been provided for this image
In [277]:
aaa_al_m2 = arima(aaa, order = c(2, 1, 3))
aaa_al_m2;
par(bg = 'white')
tsdiag(aaa_al_m2, gof = 52)
Call:
arima(x = aaa, order = c(2, 1, 3))

Coefficients:
         ar1      ar2      ma1     ma2     ma3
      1.4382  -0.6681  -1.0606  0.2126  0.2229
s.e.  0.1677   0.1017   0.1655  0.0600  0.0358

sigma^2 estimated as 0.007976:  log likelihood = 2457.81,  aic = -4905.62
No description has been provided for this image
In [288]:
aaa_al_m2 = arima(
    aaa, order = c(2, 1, 3),
    seasonal = list(order = c(1, 0, 0), period = 26)
)
aaa_al_m2;
par(bg = 'white')
tsdiag(aaa_al_m2, gof = 52)
Call:
arima(x = aaa, order = c(2, 1, 3), seasonal = list(order = c(1, 0, 0), period = 26))

Coefficients:
         ar1      ar2      ma1     ma2     ma3    sar1
      1.4202  -0.6566  -1.0427  0.2074  0.2218  0.0023
s.e.  0.1676   0.1082   0.1659  0.0632  0.0353  0.0206

sigma^2 estimated as 0.007976:  log likelihood = 2457.81,  aic = -4903.62
No description has been provided for this image

Forecast

In [282]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(9,1,0)
fixed = NULL
seasonal = NULL
aaa_fc_res = plot_arima_forecast_fig(
    da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', 
    include.mean=TRUE, transform.pars=NULL,
    main="Forecasts from ARIMA(9,1,0) for Aaa weekly bond yields", 
    xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
[1] 2440
2008.423 ; 2009.442
No description has been provided for this image
In [283]:
aaa_fc_tb = comb_forecast_res(aaa_fc_res, aaa_ts, eotr, freq)
aaa_fc_tb
Forecast method: ARIMA(9,1,0)

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
         ar1      ar2     ar3      ar4     ar5      ar6      ar7     ar8
      0.3998  -0.0872  0.1069  -0.0305  0.0836  -0.0631  -0.0200  0.0062
s.e.  0.0202   0.0218  0.0219   0.0220  0.0219   0.0220   0.0219  0.0218
          ar9
      -0.0366
s.e.   0.0203

sigma^2 estimated as 0.007577:  log likelihood = 2494.55,  aic = -4971.1

Error measures:
                       ME       RMSE        MAE         MPE      MAPE
Training set 0.0004752123 0.08702582 0.05640702 0.006786893 0.6917595
                   MASE          ACF1
Training set 0.07420383 -0.0007418544

Forecasts:
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2008.942       6.207367 6.095816 6.318917 6.036764 6.377969
2008.962       6.245664 6.053766 6.437563 5.952180 6.539148
2008.981       6.283379 6.030786 6.535972 5.897072 6.669686
2009.000       6.284329 5.976738 6.591920 5.813909 6.754749
2009.019       6.289618 5.932881 6.646355 5.744036 6.835200
2009.038       6.286249 5.881784 6.690713 5.667674 6.904823
2009.058       6.267942 5.820319 6.715565 5.583362 6.952523
2009.077       6.259608 5.774415 6.744802 5.517569 7.001648
2009.096       6.248974 5.729201 6.768747 5.454050 7.043898
2009.115       6.240264 5.689523 6.791004 5.397979 7.082548
2009.135       6.235576 5.656606 6.814546 5.350118 7.121034
2009.154       6.230780 5.625430 6.836131 5.304976 7.156585
2009.173       6.229189 5.599089 6.859289 5.265534 7.192844
2009.192       6.228524 5.574816 6.882233 5.228764 7.228285
2009.212       6.228147 5.551861 6.904432 5.193857 7.262436
2009.231       6.229018 5.531003 6.927033 5.161497 7.296540
2009.250       6.229685 5.510556 6.948814 5.129872 7.329499
2009.269       6.230454 5.490795 6.970113 5.099243 7.361665
2009.288       6.231238 5.471554 6.990922 5.069402 7.393074
2009.308       6.231713 5.452463 7.010963 5.039953 7.423473
2009.327       6.232172 5.433800 7.030544 5.011167 7.453177
2009.346       6.232437 5.415350 7.049524 4.982810 7.482064
2009.365       6.232556 5.397150 7.067963 4.954913 7.510200
2009.385       6.232639 5.379291 7.085986 4.927557 7.537720
2009.404       6.232622 5.361692 7.103553 4.900649 7.564596
2009.423       6.232587 5.344419 7.120754 4.874251 7.590922
A Time Series:
  1. 6.20736655719426
  2. 6.24566440620478
  3. 6.28337907629132
  4. 6.2843292268527
  5. 6.28961777033179
  6. 6.28624851922467
  7. 6.26794223063653
  8. 6.25960830867983
  9. 6.24897423756542
  10. 6.24026376349615
  11. 6.23557584923931
  12. 6.23078046674665
  13. 6.22918899621492
  14. 6.22852446989816
  15. 6.22814662901021
  16. 6.22901817876502
  17. 6.22968504832111
  18. 6.23045371472518
  19. 6.23123800093802
  20. 6.23171307644446
  21. 6.23217201001091
  22. 6.23243689662713
  23. 6.23255633090181
  24. 6.23263853121758
  25. 6.23262240593548
  26. 6.23258655524976
A Time Series:
  1. 0.087043610126997
  2. 0.149739441826816
  3. 0.197099142521436
  4. 0.240014699961627
  5. 0.278363324822365
  6. 0.315605045837638
  7. 0.349282180622216
  8. 0.378598464546428
  9. 0.405580825770749
  10. 0.429744933496646
  11. 0.451772454422396
  12. 0.472357863116345
  13. 0.491669908919013
  14. 0.510091289940599
  15. 0.527708392613043
  16. 0.544663892810058
  17. 0.561139660129361
  18. 0.577159126860929
  19. 0.592784567417526
  20. 0.60805182027381
  21. 0.622972964543716
  22. 0.63757659436269
  23. 0.651870987707066
  24. 0.665870342243498
  25. 0.679590630547217
  26. 0.693041042292113
A Time Series:
  1. 6.47
  2. 6.32
  3. 6.42
  4. 6.37
  5. 6.37
  6. 5.99
  7. 5.75
  8. 5.31
  9. 5.35
  10. 4.95
  11. 4.72
  12. 4.73
  13. 5.04
  14. 4.89
  15. 5.1
  16. 5.23
  17. 5.29
  18. 5.21
  19. 5.25
  20. 5.31
  21. 5.4
  22. 5.49
  23. 5.62
  24. 5.51
  25. 5.41
  26. 5.47
A Time Series: 26 × 3
ForecastStd. ErrorActual
2008.9426.2073670.087043616.47
2008.9626.2456640.149739446.32
2008.9816.2833790.197099146.42
2009.0006.2843290.240014706.37
2009.0196.2896180.278363326.37
2009.0386.2862490.315605055.99
2009.0586.2679420.349282185.75
2009.0776.2596080.378598465.31
2009.0966.2489740.405580835.35
2009.1156.2402640.429744934.95
2009.1356.2355760.451772454.72
2009.1546.2307800.472357864.73
2009.1736.2291890.491669915.04
2009.1926.2285240.510091294.89
2009.2126.2281470.527708395.10
2009.2316.2290180.544663895.23
2009.2506.2296850.561139665.29
2009.2696.2304540.577159135.21
2009.2886.2312380.592784575.25
2009.3086.2317130.608051825.31
2009.3276.2321720.622972965.40
2009.3466.2324370.637576595.49
2009.3656.2325560.651870995.62
2009.3856.2326390.665870345.51
2009.4046.2326220.679590635.41
2009.4236.2325870.693041045.47
In [293]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(2,1,3)
fixed = NULL
seasonal = NULL
aaa_fc_res_1 = plot_arima_forecast_fig(
    da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', 
    include.mean=TRUE, transform.pars=NULL,
    main="Forecasts from ARIMA(2,1,3) for Aaa weekly bond yields", 
    xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
[1] 2440
2008.423 ; 2009.442
No description has been provided for this image
In [285]:
aaa_fc_tb_1 = comb_forecast_res(aaa_fc_res_1, aaa_ts, eotr, freq)
aaa_fc_tb_1
Forecast method: ARIMA(2,1,3)

Model Information:

Call:
arima(x = tr_da_ts, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = fixed, method = method)

Coefficients:
         ar1     ar2     ma1      ma2      ma3
      0.0928  0.3274  0.3084  -0.2907  -0.0277
s.e.  0.1809  0.1336  0.1826   0.1030   0.0537

sigma^2 estimated as 0.007631:  log likelihood = 2485.79,  aic = -4961.58

Error measures:
                       ME       RMSE        MAE         MPE      MAPE
Training set 0.0004452245 0.08733975 0.05630836 0.006843011 0.6904256
                   MASE         ACF1
Training set 0.07407404 2.510243e-05

Forecasts:
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2008.942       6.203303 6.091350 6.315257 6.032086 6.374521
2008.962       6.222434 6.029719 6.415148 5.927702 6.517165
2008.981       6.247466 5.993676 6.501256 5.859328 6.635605
2009.000       6.256052 5.946342 6.565763 5.782391 6.729714
2009.019       6.265046 5.906135 6.623956 5.716140 6.813951
2009.038       6.268691 5.864548 6.672835 5.650608 6.886775
2009.058       6.271974 5.826473 6.717475 5.590639 6.953309
2009.077       6.273473 5.789517 6.757428 5.533327 7.013618
2009.096       6.274687 5.754862 6.794512 5.479683 7.069690
2009.115       6.275290 5.721714 6.828866 5.428668 7.121912
2009.135       6.275743 5.690267 6.861220 5.380334 7.171153
2009.154       6.275983 5.660188 6.891778 5.334206 7.217760
2009.173       6.276154 5.631432 6.920875 5.290137 7.262170
2009.192       6.276248 5.603819 6.948677 5.247857 7.304639
2009.212       6.276313 5.577262 6.975363 5.207207 7.345418
2009.231       6.276350 5.551646 7.001053 5.168012 7.384687
2009.250       6.276374 5.526891 7.025857 5.130139 7.422609
2009.269       6.276388 5.502916 7.049861 5.093465 7.459312
2009.288       6.276398 5.479657 7.073139 5.057888 7.494908
2009.308       6.276403 5.457053 7.095754 5.023315 7.529491
2009.327       6.276407 5.435054 7.117760 4.989669 7.563145
2009.346       6.276409 5.413614 7.139204 4.956878 7.595940
2009.365       6.276411 5.392694 7.160127 4.924882 7.627939
2009.385       6.276411 5.372256 7.180567 4.893625 7.659197
2009.404       6.276412 5.352270 7.200553 4.863059 7.689765
2009.423       6.276412 5.332707 7.220117 4.833140 7.719684
A Time Series:
  1. 6.20330330269006
  2. 6.22243353443424
  3. 6.24746607245773
  4. 6.25605234307303
  5. 6.26504558322429
  6. 6.26869137771544
  7. 6.27197435722184
  8. 6.27347270037949
  9. 6.27468668481252
  10. 6.27528992475406
  11. 6.2757433959077
  12. 6.27598298982167
  13. 6.276153702089
  14. 6.27624799191435
  15. 6.27631263740499
  16. 6.27634950891711
  17. 6.27637409708012
  18. 6.27638845136805
  19. 6.27639783419081
  20. 6.27640340481769
  21. 6.27640699392653
  22. 6.27640915093988
  23. 6.27641052626771
  24. 6.27641136015153
  25. 6.27641188784976
  26. 6.27641220985211
A Time Series:
  1. 0.0873575969045483
  2. 0.150376034883939
  3. 0.198033479167829
  4. 0.241668421990196
  5. 0.280059077131007
  6. 0.315354557154003
  7. 0.347626222826601
  8. 0.377632255950079
  9. 0.405621476362031
  10. 0.431957844624765
  11. 0.45684997407101
  12. 0.480507067616555
  13. 0.503078925929186
  14. 0.524698835464767
  15. 0.545471993235578
  16. 0.565488940301752
  17. 0.584824634122161
  18. 0.603543457479824
  19. 0.621700203437587
  20. 0.63934233836289
  21. 0.656511023822344
  22. 0.673242320749988
  23. 0.68956795117004
  24. 0.705516016121785
  25. 0.721111524730885
  26. 0.736376854489315
A Time Series:
  1. 6.47
  2. 6.32
  3. 6.42
  4. 6.37
  5. 6.37
  6. 5.99
  7. 5.75
  8. 5.31
  9. 5.35
  10. 4.95
  11. 4.72
  12. 4.73
  13. 5.04
  14. 4.89
  15. 5.1
  16. 5.23
  17. 5.29
  18. 5.21
  19. 5.25
  20. 5.31
  21. 5.4
  22. 5.49
  23. 5.62
  24. 5.51
  25. 5.41
  26. 5.47
A Time Series: 26 × 3
ForecastStd. ErrorActual
2008.9426.2033030.08735766.47
2008.9626.2224340.15037606.32
2008.9816.2474660.19803356.42
2009.0006.2560520.24166846.37
2009.0196.2650460.28005916.37
2009.0386.2686910.31535465.99
2009.0586.2719740.34762625.75
2009.0776.2734730.37763235.31
2009.0966.2746870.40562155.35
2009.1156.2752900.43195784.95
2009.1356.2757430.45685004.72
2009.1546.2759830.48050714.73
2009.1736.2761540.50307895.04
2009.1926.2762480.52469884.89
2009.2126.2763130.54547205.10
2009.2316.2763500.56548895.23
2009.2506.2763740.58482465.29
2009.2696.2763880.60354355.21
2009.2886.2763980.62170025.25
2009.3086.2764030.63934235.31
2009.3276.2764070.65651105.40
2009.3466.2764090.67324235.49
2009.3656.2764110.68956805.62
2009.3856.2764110.70551605.51
2009.4046.2764120.72111155.41
2009.4236.2764120.73637695.47
In [298]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(2,1,3)
fixed = NULL
seasonal = list(order = c(0, 0, 1), period = 13)
aaa_fc_res_2 = plot_arima_forecast_fig(
    da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq, 
    order=order, seasonal=seasonal, fixed=fixed, method='ML', 
    include.mean=TRUE, transform.pars=NULL,
    main="Forecasts from ARIMA(2,1,3)(0,0,1)[13] for Aaa weekly bond yields", 
    xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
[1] 2440
2008.423 ; 2009.442
No description has been provided for this image
In [ ]:

2-12

In [299]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da2 = read.table("../AFTS_sol/data/w-Baa.txt", header = FALSE)
da1[1:5,]; da2[1:5,];
A data.frame: 5 × 4
V1V2V3V4
<int><int><int><dbl>
119621 54.43
219621124.42
319621194.42
419621264.41
519622 24.42
A data.frame: 5 × 4
V1V2V3V4
<int><int><int><dbl>
119621 55.11
219621125.09
319621195.08
419621265.08
519622 25.07
In [302]:
freq = 52
start = c(1962, 1)
aaa = da1[,4]
aaa_ts = ts(aaa, frequency = freq, start = start)
bbb = da2[,4]
bbb_ts = ts(bbb, frequency = freq, start = start)
aaac = diff(aaa)
bbbc = diff(bbb)
aaac_ts = ts(aaac, frequency = freq, start = start)
bbbc_ts = ts(bbbc, frequency = freq, start = start)
In [316]:
apply(cbind(aaac, bbbc), 2, basicStats)
$aaac
A data.frame: 16 × 1
X..newX..i
<dbl>
nobs2466.000000
NAs 0.000000
Minimum -0.700000
Maximum 0.550000
1. Quartile -0.040000
3. Quartile 0.040000
Mean 0.000422
Median 0.000000
Sum 1.040000
SE Mean 0.001946
LCL Mean -0.003394
UCL Mean 0.004238
Variance 0.009339
Stdev 0.096638
Skewness -0.654390
Kurtosis 8.061043
$bbbc
A data.frame: 16 × 1
X..newX..i
<dbl>
nobs2466.000000
NAs 0.000000
Minimum -0.680000
Maximum 0.840000
1. Quartile -0.040000
3. Quartile 0.040000
Mean 0.001407
Median 0.000000
Sum 3.470000
SE Mean 0.001714
LCL Mean -0.001954
UCL Mean 0.004768
Variance 0.007244
Stdev 0.085113
Skewness 0.133853
Kurtosis 9.247167

Exploration

In [301]:
par(bg = 'white')
plot(aaa_ts, main = "Weekly Aaa or Bbb Bond Yield", xlab = "Year", ylab = "Yield")
lines(bbb_ts, lty = 2)
legend(
    "topleft",
    legend = c("Aaa", "Bbb"),
    lty = c(1, 2)
)
No description has been provided for this image
In [304]:
par(bg = 'white')
plot(aaac_ts, main = "Weekly Aaa Bond Yield", xlab = "Year", ylab = "Yield")
plot(bbbc_ts, main = "Weekly Bbb Bond Yield", xlab = "Year", ylab = "Yield")
No description has been provided for this image
No description has been provided for this image
In [305]:
par(bg = "white")
plot(x = aaa, y = bbb, main = "Aaa vs Bbb", xlab = "Aaa Yield", ylab = "Bbb Yield")
plot(x = aaac, y = bbbc, main = "Aaa Change vs Bbb Change", xlab = "diff(Aaa Yield)", ylab = "diff(Bbb Yield)")
No description has been provided for this image
No description has been provided for this image

Fitting

  • Strong serial correlation in ACF of residuals, indicating that there is a unit-root process.
  • Also conducted ADF unit-root tests. The results show that null-hypothesis cannot be rejected.
    • Null-hypothesis is that the time series has unit-root.
  • Box-Ljung test shows that numm-hypothesis can be rejected.
    • Null-hypothesis is that the time series at certain lag has no autocorrelation.
  • After taking one diff, the autocorrelation is largely gone. But Box-Ljung test is still significant. We can model the residuals as a time series.

Fitting the bond yields

In [306]:
bond_m1 = lm(aaa~bbb)
summary(bond_m1)
Call:
lm(formula = aaa ~ bbb)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.46094 -0.14265  0.06717  0.20420  0.68515 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.030569   0.023025   1.328    0.184    
bbb         0.881591   0.002488 354.350   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3357 on 2465 degrees of freedom
Multiple R-squared:  0.9807,	Adjusted R-squared:  0.9807 
F-statistic: 1.256e+05 on 1 and 2465 DF,  p-value: < 2.2e-16
In [308]:
par(bg = 'white')
resi_ts = ts(bond_m1$residuals, frequency = freq, start = start)
plot(resi_ts, type = 'l', main = "Simple Linear Model", xlab = "Year")
lag.max = 36
par(mfrow = c(2,1))
pacf(bond_m1$residuals, lag.max = lag.max, main = "Residuals PACF Plot")
acf(bond_m1$residuals, lag.max = lag.max, main = "Residuals ACF Plot")
No description has been provided for this image
No description has been provided for this image
In [309]:
adf_test_res_a = adfTest(aaa, lags = 12, type = c("ct"))
adf_test_res_b = adfTest(bbb, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_a; adf_test_res_b
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -1.592
  P VALUE:
    0.7511 

Description:
 Thu Feb 29 12:01:44 2024 by user: 
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -1.6934
  P VALUE:
    0.7081 

Description:
 Thu Feb 29 12:01:44 2024 by user: 
In [310]:
Box.test(bond_m1$residuals, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  bond_m1$residuals
X-squared = 21007, df = 12, p-value < 2.2e-16

Fitting the bond yield changes

In [311]:
bond_m2 = lm(aaac~-1+bbbc) # No intercept
summary(bond_m2)
Call:
lm(formula = aaac ~ -1 + bbbc)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.44475 -0.01838  0.00000  0.01946  0.36192 

Coefficients:
     Estimate Std. Error t value Pr(>|t|)    
bbbc  0.94613    0.01264   74.87   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05341 on 2465 degrees of freedom
Multiple R-squared:  0.6946,	Adjusted R-squared:  0.6944 
F-statistic:  5605 on 1 and 2465 DF,  p-value: < 2.2e-16
In [312]:
par(bg = 'white')
resi_ts = ts(bond_m2$residuals, frequency = freq, start = start)
plot(resi_ts, type = 'l', main = "Simple Linear Model", xlab = "Year")
lag.max = 36
par(mfrow = c(2,1))
pacf(bond_m2$residuals, lag.max = lag.max, main = "Residuals PACF Plot")
acf(bond_m2$residuals, lag.max = lag.max, main = "Residuals ACF Plot")
No description has been provided for this image
No description has been provided for this image
In [314]:
adf_test_res_ac = adfTest(aaac, lags = 12, type = c("ct"))
adf_test_res_bc = adfTest(bbbc, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_ac; adf_test_res_bc
Warning message in adfTest(aaac, lags = 12, type = c("ct")):
“p-value smaller than printed p-value”
Warning message in adfTest(bbbc, lags = 12, type = c("ct")):
“p-value smaller than printed p-value”
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -14.495
  P VALUE:
    0.01 

Description:
 Thu Feb 29 12:09:09 2024 by user: 
Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -13.0501
  P VALUE:
    0.01 

Description:
 Thu Feb 29 12:09:09 2024 by user: 
In [313]:
Box.test(bond_m2$residuals, lag = 12, type = 'Ljung')
	Box-Ljung test

data:  bond_m2$residuals
X-squared = 189.69, df = 12, p-value < 2.2e-16

Fitting a time series residuals

In [319]:
bond_m3 = arima(aaac, order = c(0,0,1), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(0, 0, 1), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ma1    xreg
      0.2335  0.9436
s.e.  0.0185  0.0132

sigma^2 estimated as 0.00269:  log likelihood = 3797.81,  aic = -7591.62

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.628555479175029
No description has been provided for this image
In [320]:
bond_m3 = arima(aaac, order = c(1,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(1, 0, 0), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1    xreg
      0.2394  0.9472
s.e.  0.0195  0.0133

sigma^2 estimated as 0.002688:  log likelihood = 3798.99,  aic = -7593.98

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.628910358713633
No description has been provided for this image
In [323]:
bond_m3 = arima(aaac, order = c(1,0,1), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(1, 0, 1), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1     ma1    xreg
      0.1476  0.0980  0.9455
s.e.  0.0710  0.0706  0.0133

sigma^2 estimated as 0.002686:  log likelihood = 3799.9,  aic = -7593.81

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.629186078699108
No description has been provided for this image
In [ ]:

In [324]:
cat("order =", ar(bond_m2$residuals, method = 'mle')$order, "\n")
order = 12 
In [325]:
bond_m3 = arima(aaac, order = c(12,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(12, 0, 0), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1      ar2      ar3     ar4     ar5      ar6      ar7      ar8
      0.2430  -0.0208  -0.0432  0.0632  0.0233  -0.0375  -0.0388  -0.0215
s.e.  0.0201   0.0208   0.0207  0.0208  0.0208   0.0208   0.0209   0.0208
          ar9     ar10     ar11    ar12    xreg
      -0.0129  -0.0235  -0.0451  0.0294  0.9409
s.e.   0.0208   0.0207   0.0207  0.0202  0.0133

sigma^2 estimated as 0.002646:  log likelihood = 3818.08,  aic = -7610.16

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.634626474385314
No description has been provided for this image
In [326]:
bond_m3 = arima(aaac, order = c(2,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(2, 0, 0), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1      ar2    xreg
      0.2468  -0.0310  0.9453
s.e.  0.0201   0.0202  0.0133

sigma^2 estimated as 0.002685:  log likelihood = 3800.16,  aic = -7594.32

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.629263915621719
No description has been provided for this image
In [327]:
bond_m3 = arima(aaac, order = c(3,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(3, 0, 0), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1      ar2      ar3    xreg
      0.2459  -0.0245  -0.0257  0.9463
s.e.  0.0201   0.0208   0.0202  0.0133

sigma^2 estimated as 0.002684:  log likelihood = 3800.97,  aic = -7593.94

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.629507134709869
No description has been provided for this image
In [338]:
bond_m3 = arima(aaac, order = c(4,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(4, 0, 0), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1      ar2      ar3     ar4    xreg
      0.2477  -0.0226  -0.0430  0.0698  0.9475
s.e.  0.0201   0.0208   0.0207  0.0201  0.0132

sigma^2 estimated as 0.002671:  log likelihood = 3806.98,  aic = -7603.96

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.631311205650992
No description has been provided for this image
In [339]:
bond_m3 = arima(aaac, order = c(5,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(5, 0, 0), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1      ar2      ar3     ar4     ar5    xreg
      0.2469  -0.0222  -0.0427  0.0672  0.0107  0.9475
s.e.  0.0201   0.0208   0.0207  0.0207  0.0202  0.0133

sigma^2 estimated as 0.00267:  log likelihood = 3807.12,  aic = -7602.24

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.631353074608109
No description has been provided for this image
In [340]:
bond_m3 = arima(aaac, order = c(8,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(8, 0, 0), xreg = bbbc, include.mean = FALSE)

Coefficients:
         ar1      ar2      ar3     ar4     ar5      ar6      ar7      ar8
      0.2443  -0.0194  -0.0410  0.0652  0.0215  -0.0408  -0.0385  -0.0226
s.e.  0.0201   0.0208   0.0208  0.0207  0.0207   0.0207   0.0208   0.0202
        xreg
      0.9424
s.e.  0.0133

sigma^2 estimated as 0.002657:  log likelihood = 3813.3,  aic = -7608.6

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.633201591685739
No description has been provided for this image
In [335]:
bond_m3 = arima(aaac, order = c(2,0,2), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(2, 0, 2), xreg = bbbc, include.mean = FALSE)

Coefficients:
          ar1      ar2     ma1     ma2    xreg
      -0.1971  -0.2378  0.4430  0.3195  0.9469
s.e.   0.2242   0.1148  0.2255  0.0810  0.0132

sigma^2 estimated as 0.002675:  log likelihood = 3805.1,  aic = -7600.21

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.630749356776852
No description has been provided for this image
In [337]:
bond_m3 = arima(
    aaac, order = c(2,0,2), 
    fixed = c(0, NA, NA, NA, NA),
    xreg = bbbc, include.mean = FALSE
)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Warning message in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg, :
“some AR parameters were fixed: setting transform.pars = FALSE”
Warning message in trainingaccuracy(object, test, d, D):
“test elements must be within sample”
Call:
arima(x = aaac, order = c(2, 0, 2), xreg = bbbc, include.mean = FALSE, fixed = c(0, 
    NA, NA, NA, NA))

Coefficients:
      ar1      ar2     ma1     ma2    xreg
        0  -0.2829  0.2464  0.3113  0.9464
s.e.    0   0.0982  0.0205  0.0925  0.0131

sigma^2 estimated as 0.002676:  log likelihood = 3804.53,  aic = -7601.07

Training set error measures:
              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN
0.63057831151478
No description has been provided for this image

Forecast

In [332]:
npts = 26
eotr = length(aaac_ts)-npts
h = npts
freq = 52
xreg = as.matrix(bbbc)
aaac_fc_res = plot_auto_arima_forecast_fig(
    da_ts=aaac_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
    xreg=xreg,
    main=NULL, # "Forecasts from ARIMA(2,0,2)\nfor CRSP Equal-Weighted Index",
    xlab="Year", ylab="Aaa Bond Yield (using Baa Bond Yield)", # , ylim=c(5, 11)
    d = 0,
    D = 0,
    max.p = 3,
    max.q = 4,
    max.P = 0,
    max.Q = 0,
    max.order = 5,
    seasonal = FALSE,
    method = "ML",
    allowmean = FALSE,
    stepwise = FALSE,
    parallel = TRUE,
    num.cores = 2
)
[1] 2440
2008.404 ; 2009.423
No description has been provided for this image
In [333]:
aaac_fc_tb = comb_forecast_res(aaac_fc_res, aaac_ts, eotr, freq)
aaac_fc_tb
Forecast method: Regression with ARIMA(2,0,3) errors

Model Information:
Series: tr_da_ts 
Regression with ARIMA(2,0,3) errors 

Coefficients:
          ar1      ar2     ma1     ma2     ma3    xreg
      -0.0073  -0.9447  0.2120  0.9378  0.1590  0.9570
s.e.   0.0154   0.0361  0.0253  0.0396  0.0236  0.0131

sigma^2 = 0.002507:  log likelihood = 3846.97
AIC=-7679.95   AICc=-7679.9   BIC=-7639.35

Error measures:
                        ME       RMSE        MAE MPE MAPE     MASE        ACF1
Training set -0.0004666278 0.05000607 0.03054997 NaN  Inf 0.344656 0.006556348

Forecasts:
         Point Forecast        Lo 80        Hi 80        Lo 95       Hi 95
2008.923    0.768124723  0.703960422  0.832289023  0.669993894  0.86625555
2008.942    0.196099159  0.130605234  0.261593084  0.095934844  0.29626347
2008.962    0.195199605  0.129703479  0.260695730  0.095031925  0.29536728
2008.981   -0.157585923 -0.223118946 -0.092052899 -0.257810033 -0.05736181
2009.000   -0.070542178 -0.136077291 -0.005007065 -0.170769484  0.02968513
2009.019   -0.110605586 -0.176173488 -0.045037684 -0.210883039 -0.01032813
2009.038   -0.034959451 -0.100529333  0.030610431 -0.135239932  0.06532103
2009.058   -0.310281374 -0.375880393 -0.244682355 -0.410606416 -0.20995633
2009.077   -0.060531180 -0.126132072  0.005069712 -0.160859086  0.03979673
2009.096   -0.398127384 -0.463754165 -0.332500603 -0.498494884 -0.29775988
2009.115   -0.198069283 -0.263697832 -0.132440733 -0.298439487 -0.09769908
2009.135   -0.022779233 -0.088430784  0.042872318 -0.123184616  0.07762615
2009.154    0.150404876  0.084751658  0.216058094  0.049996943  0.25081281
2009.173   -0.245372922 -0.311046575 -0.179699268 -0.345812107 -0.14493374
2009.192    0.174812481  0.109137258  0.240487703  0.074370895  0.27525407
2009.212    0.063708105 -0.001985271  0.129401481 -0.036761243  0.16417745
2009.231   -0.040662205 -0.106357056  0.025032645 -0.141133809  0.05980940
2009.250   -0.159576203 -0.225287178 -0.093865228 -0.260072467 -0.05907994
2009.269    0.002226129 -0.063486230  0.067938489 -0.098272252  0.10272451
2009.288    0.111880532  0.046153851  0.177607213  0.011360248  0.21240082
2009.308    0.093622583  0.027894605  0.159350562 -0.006899685  0.19414485
2009.327    0.165512032  0.099771335  0.231252729  0.064970312  0.26605375
2009.346    0.097649441  0.031907530  0.163391353 -0.002894136  0.19819302
2009.365    0.026037185 -0.039716021  0.091790390 -0.074523665  0.12659803
2009.385   -0.059240699 -0.124995040  0.006513642 -0.159803286  0.04132189
2009.404    0.107813680  0.042049311  0.173578049  0.007235757  0.20839160
A Time Series:
  1. 0.768124722892773
  2. 0.196099158934993
  3. 0.195199604711239
  4. -0.157585922811906
  5. -0.0705421780087563
  6. -0.110605585759758
  7. -0.0349594511087951
  8. -0.31028137413802
  9. -0.060531180049755
  10. -0.398127383690043
  11. -0.198069282747511
  12. -0.0227792332303308
  13. 0.15040487594322
  14. -0.245372921535585
  15. 0.17481248068782
  16. 0.0637081051064078
  17. -0.0406622054036674
  18. -0.159576202991684
  19. 0.00222612936507818
  20. 0.111880531657912
  21. 0.093622583315119
  22. 0.165512031566735
  23. 0.0976494413530556
  24. 0.0260371848928738
  25. -0.0592406989075704
  26. 0.107813679982687
A Time Series:
  1. 0.0500676696689963
  2. 0.0511051812392695
  3. 0.0511068983526645
  4. 0.0511356898710673
  5. 0.0511373202138411
  6. 0.0511629059237635
  7. 0.0511644508783508
  8. 0.0511871863791592
  9. 0.0511886477400432
  10. 0.0512088491938431
  11. 0.051210229095681
  12. 0.0512281776296255
  13. 0.0512294784825296
  14. 0.0512454240904818
  15. 0.0512466485221148
  16. 0.0512608135248649
  17. 0.0512619643289661
  18. 0.0512745464447021
  19. 0.0512756265356275
  20. 0.0512868016059121
  21. 0.051287813978761
  22. 0.0512977383897053
  23. 0.0512986860855978
  24. 0.0513074989191092
  25. 0.0513083849951865
  26. 0.0513162099409904
A Time Series:
  1. 0.35
  2. -0.149999999999999
  3. 0.0999999999999996
  4. -0.0499999999999998
  5. 0
  6. -0.38
  7. -0.24
  8. -0.44
  9. 0.04
  10. -0.399999999999999
  11. -0.23
  12. 0.0100000000000007
  13. 0.31
  14. -0.15
  15. 0.21
  16. 0.130000000000001
  17. 0.0599999999999996
  18. -0.0800000000000001
  19. 0.04
  20. 0.0599999999999996
  21. 0.0900000000000007
  22. 0.0899999999999999
  23. 0.13
  24. -0.11
  25. -0.0999999999999996
  26. 0.0599999999999996
A Time Series: 26 × 3
ForecastStd. ErrorActual
2008.923 0.7681247230.05006767 0.35
2008.942 0.1960991590.05110518-0.15
2008.962 0.1951996050.05110690 0.10
2008.981-0.1575859230.05113569-0.05
2009.000-0.0705421780.05113732 0.00
2009.019-0.1106055860.05116291-0.38
2009.038-0.0349594510.05116445-0.24
2009.058-0.3102813740.05118719-0.44
2009.077-0.0605311800.05118865 0.04
2009.096-0.3981273840.05120885-0.40
2009.115-0.1980692830.05121023-0.23
2009.135-0.0227792330.05122818 0.01
2009.154 0.1504048760.05122948 0.31
2009.173-0.2453729220.05124542-0.15
2009.192 0.1748124810.05124665 0.21
2009.212 0.0637081050.05126081 0.13
2009.231-0.0406622050.05126196 0.06
2009.250-0.1595762030.05127455-0.08
2009.269 0.0022261290.05127563 0.04
2009.288 0.1118805320.05128680 0.06
2009.308 0.0936225830.05128781 0.09
2009.327 0.1655120320.05129774 0.09
2009.346 0.0976494410.05129869 0.13
2009.365 0.0260371850.05130750-0.11
2009.385-0.0592406990.05130838-0.10
2009.404 0.1078136800.05131621 0.06
In [342]:
par(bg = 'white')
tsdiag(aaac_fc_res$model, gof.lag = 36)
No description has been provided for this image
In [ ]:

2-13

In [ ]:

In [ ]:

In [ ]:

In [ ]:

2-14

In [ ]:

In [ ]:

In [ ]:

In [ ]:

2-15

In [ ]:

In [ ]:

In [ ]: